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ABSTRACT 


Some  new  theory  and  applications  of  certain  bivariate 
non-normal  distributions  are  presented.  In  particular,  new 
bivariate  negative  binomial  and  gamma  distributions  are  dis- 
cussed and  an  existing  bivariate  exponential  distribution  is 
applied  to  single  stage  and  tandem  queueing  systems  (both 
single  server)  which  have  particular  kinds  of  correlated 
structure. 

One  new  bivariate  negative  binomial  distribution  is  de- 
rived by  convoluting  an  existing  bivariate  geometric  distri- 
bution; the  probability  function  has  six  parameters  and  ad- 
mits of  positive  or  negative  correlations  and  linear  or  non- 
linear regressions.  Given  are  the  moments  to  order  two  and 
for  special  cases,  the  regression  function,  a recursive  form- 
ula for  the  probabilities,  a method  of  moments  parameter  es- 
timation technique,  the  likelihood  equations,  the  differential 
difference  equations  and  for  maximum  likelihood  parameter  es- 
timates, a necessary  relationship  for  the  parameters.  Certain 
results  are  extended  to  a dual  bivariate  gamma  distribution. 
Another  bivariate  negative  binomial  distribution,  which  has 
four  parameters,  results  by  reducing  a particular  trivariate 
negative  binomial  distribution  with  independent  marginals; 
only  positive  correlations  and  linear  regressions  are  possible 
here.  Both  bivariate  negative  binomial  distributions  are 


x 


If 
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fitted  to  data  and  their  special  features  illustrated. 

Applications  of  bivariate  distributions  to  certain  air- 
craft logistical  problems  are  investigated.  Primarily,  a 
bivariate  negative  binomial  distribution  is  fitted  to  spare 
parts  demand  data  in  two  periods  and  to  monthly  abort  data 
on  either  side  of  a large  scale  maintenance  event  and  it  is 
shown  how  the  associated  sample  distributions  can  be  useful 
in  parts  inventory  control  and  in  investigating  the  effect 
of  maintenance  on  an  aircraft's  performance. 

A bivariate  exponential  distribution  is  applied  to  tan- 
dem queues  to  study  the  effect  of  correlated  exponential  ser- 
vice times  and  to  single  stage  queues  to  study  the  effect  of 
correlation  between  a customer's  service  time  and  the  inter- 
arrival interval  separating  himself  and  his  predecessor.  Ar- 
rivals to  both  systems  are  according  to  a Poisson  process. 
Simulation  is  used  to  show  that  the  mean  waiting  time  is  quite 
sensitive  to  departures  from  the  traditional  assumptions  of 
mutually  independent  service  times  for  tandem  queues  and  inde- 
pendence of  service  times  and  interarrival  intervals  for  single 
stage  queues,  especially  at  higher  utilizations.  For  the  cases 
of  infinite  interstage  storage  between  two-stage  tandem  queues 
and  infinite  storage  before  a single  stage  queue, system  per- 
formance is  increased  by  positive  correlation  and  impaired  by 

% 

negative  correlation.  For  two-stage  queues  this  change  is  re- 
versed for  zero  interstage  storage  and  depends  on  the  value  of 
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the  utilization  rate  for  the  case  where  interstage  storage 
equals  unity.  By  using  spectral  analysis  techniques  and  a 
nonparametric  test  applied  to  sample  power  spectra  associ- 
ated with  certain  simulated  waiting  times  the  effects  are 
shown  to  be  statistically  significant.  For  correlation 
equal  unity  and  infinite  interstage  storage  results  are 
given  for  two  through  twenty-five  stages  in  series. 


PART  1 


INTRODUCTION 

The  objectives  of  this  research  are  (1)  to  develop 
two  new  bivariate  negative  binomial  probability  functions, 

(2)  to  derive,  where  applicable,  corresponding  properties  for 
a dual  bivariate  gamma  distribution,  and  (3)  to  show  that  bi- 
variate approaches  to  data  analysis  and  mathematical  modeling 
can  provide,  in  some  cases,  more  meaningful  and  representative 
results  than  traditional  approaches.  The  first  bivariate  nega- 
tive binomial  (bnb)  distribution  is  derived,  via  convolution, 
from  an  existing  bivariate  geometric  distribution  and  the  second 
one  is  developed,  via  reduction,  from  a certain  independent  tri- 
variate negative  binomial  distribution.  Certain  properties  re- 
lated to  infinite  divisibility  and  parameter  estimation  are 
shown  to  be  directly  applicable  to  an  existing  bivariate  gamma 
distribution. 

In  analyzing  bivariate  data  from  self-pairing  type 
studies  usually  the  data  are  transformed  to  obtain  a univariate 
random  variable;  sometimes  this  technique  is  acceptable  partic- 
ularly if  the  data  are  approximately  normal  but  often  a bivariate 
random  variable,  say,  representing  count  data,  cannot  be  ade- 
quately treated  in  this  way  and  so  an  alternate  approach  is 
deemed  necessary.  Most  often  in  these  cases  a first  step  is 
to  find  a reasonable  bivariate  probability  function  to  repre- 
sent the  data  and  for  this  purpose  several  probability  functions 
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are  usually  compared.  That  these  new  bnb  distributions  should 
be  useful  in  practical  applications  is  suggested  by  showing, 
for  two  bivariate  samples  from  the  literature,  how  certain  pro- 
perties of  these  distributions  better  represent  the  data.  Addi- 
tionally, we  introduce  some  new  bivariate  data  sets  related  to 
aircraft  operations  and  maintenance  and  show  how  bivariate  ap- 
proaches can  be  useful  in  analyzing  certain  problems  dealing 
with  these  data. 

Another  new  bivariate  approach  is  related  to  corre- 
lated queueing  systems.  For  instance,  single  server  tandem 
queues  traditionally  have  been  modeled  by  assuming  that  cus- 
tomer service  times  at  the  individual  servers  are  independent; 
sometimes  this  is  a reasonable  assumption  but  in  certain  impor- 
tant applications,  for  example,  production  lines,  such  is  often 
not  the  case  and  a more  realistic  model  is  desired.  We  show 
for  two  servers  in  series  the  effect  of  correlated  exponential 
service  times  by  assuming  that  a customer's  service  times  at 
the  two  servers  are  given  by  a bivariate  exponential  distribu- 
tion. Also  we  provide  a similar  analysis  for  single  server, 
single  stage  queues  with  correlated  interarrival  and  service 
processes . 

Part  2 describes  the  new  bnb  distributions  with  cer- 
tain associated  results  for  a bivariate  gamma  distribution  and 
also  shows  how  the  bnb  distributions  fit  some  empirical  data. 

We  introduce  in  Part  3 some  new  discrete  bivariate  data  sets 


related  to  aircraft  operations  and  maintenance  and  show  how 
bnb  distributions  can  be  useful  in  analyzing  certain  problems 
dealing  with  these  data.  Also  we  describe  how  bivariate  gamma 
distributions  could  be  useful  in  similar  problems  associated 
with  continuous  data.  Part  4 shows  the  results  on  correlated 
queueing  systems.  More  introductory  and  historical  remarks 
are  included  in  these  parts. 


PART  2 

ON  BIVARIATE  NEGATIVE  BINOMIAL 
(AND  GAMMA)  DISTRIBUTIONS 


2.1  Introduction  and  Historical  Review 


In  this  part  we  develop  and  fit  to  data  two  new  bnb  dis- 
tributions and  show  how  certain  properties  relate  to  an  existing 
bivariate  gamma  distribution. 

That  a bnb  should  be  important  in  statistical  theory  and 
applications  is  suggested  by  the  wide  acceptance  of  the  uni- 
variate negative  binomial  distribution  as  a reasonable  model 
for  a broad  range  of  problems  representing  univariate  discrete 
random  variables  (see  Boswell  and  Patil  (1970)  foT  a discussion 
of  fifteen  stochastic  models  which  give  rise  to  the  univariate 
negative  binomial  distribution) . It  is  reasonable  to  suspect 
that  a bnb  distribution  would  be  useful  in  describing  bivariate 
random  variables  for  which  correlation  exist  between  the  members 
of  the  bivariate  pair  and  the  marginals  are  negative  binomial. 

A few  bnb  distributions  have  been  presented  in  the  literature; 
next  we  show  the  particular  forms  of  the  univariate  negative 
binomial  to  be  used  herein  and  then  review  these  bnb  distribu- 
tions and  some  of  their  applications. 

The  univariate  negative  binomial  with  parameters  v>0  and 
6>0  is  defined  (Johnson  and  Kotz  (1969))  as  the  distribution  of 
a random  variable  (r.v.)  X for  which 
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PrtX=xJ  " xffv-1  j i (T^e)V(lTU)X'  x®°  * 1 » 2 * . . 
and  the  characteristic  function  is 


(2.1) 


♦ (t)  - EteltX]  * [l+0(l-eu)]'v. 


it>i  -v 


(2.2) 


The  mean  and  variance  of  X are  v6  and  v0(l+0),  respectively. 
Thus,  it  is  characteristic  of  the  negative  binomial  distribu- 
tion that  the  variance  is  greater  than  the  mean.  A method  of 
moments  parameter  estimation  technique  is  described  in  Johnson 
and  Kotz.  Par  v*=l  we  have  the  geometric  distribution. 

Another  common  representation  is  to  let  0=p/(l-p),  or 
equivalently,  p=0/(l+0) , in  (2.1)  and  so 

Pr[X“xl  = xT^T-'l]  ! (1*P)V  PX»  *-0,1,2,...  (2.3) 


and 


♦to  - u.jE-a-e11)]-'’  - [ijKilr'’ 


(2.4) 


This  latter  representation  is  referred  to  as  a negative  bi- 
nomial distribution  with  parameters  v and  p.  We  use  both 
forms  throughout. 

The  negative  binomial  distribution  can  be  viewed  as  a 
compound  distribution  (Johnson  and  Kotz).  In  fact,  a mixture 
of  Poisson  distributions  such  that  the  expected  values,  A,  of 
the  Poisson  distributions  vary  according  to  a Type  III  (gamma) 
distribution  with  probability  density  function 
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f(^)  “ (-g-)V  \>0;  v>0,  9>0  (2.5) 

leads  to  2.1.  A multivariate  negative  binomial  distribution 
was  constructed  in  an  analogous  way  by  Bates  and  Neyman  (1952)  . 
We  describe  the  bivariate  case.  Suppose  we  consider  the  joint 
distribution  of  the  two  independent  random  variables  X and  Y 
where  X is  distributed  as  a Poisson  r.v.  with  expected  value  A 
and  Y is  also  a Poisson  r.v.  but  its  expected  value  is  aA , a>0 
and  constant.  If  we  assume  that  A is  distributed  according  to 
the  gamma  distribution  in  (2.5),  then  the  marginal  distribution 
of  X and  Y is 


Pr  pt  *x,Y«y]  - <1“P‘<l)V  pV".  x(y-0 ,1 ,2  , . . . (2.6) 


where  p*0/ [1+  (ot+1)  0]  , q-ap,  v>0,  0 <p  <1,0  <q  <1,  and 
0 < p+q  < 1.  The  characteristic  function  is 


♦(tx,t2)  - E[eitlX»it2Yl 
or  in  terms  of  0, 

♦(t1,t2)  - [l+0(l-eltl)  ♦ o0(l-eit2)] _v.  (2.8) 

We  have  that  the  mean  vector  is 


(2.9) 


V 

v6 

y ■ 

■ 

WY 

va6 

and  the  covariance  matrix  is 
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°XY 

°XY  °Y 


ve(i+e)  vote2 
2 

va9  va0(l+a0) 


(2.10) 


Prom  the  characteristic  function  it  is  clear  that  the  margin- 
als are  again  negative  binomial.  The  conditional  distribution 
can  be  shown  to  be 


Pr[Y|x] 


_ (x+y+v-1) ) 
y !'(x+v-l) ! 


d-q)X+V  qy, 


y-0,1,2, . . . 


(2.11) 


or  a negative  binomial  with  parameters  x+v  and  q.  Therefore, 
the  conditional  mean  or  regression  function  is 


E [Y | x]  = q (v+x) / (1-q) 


(2.12) 


and  we  note  that  the  form  is  linear.  Note  also  that  this 
probability  function  admits  of  positive  correlations  only. 

Besides  Bates  and  Neyman  in  1952  others  have  studied  the 
above  bnb  distribution  (Mardia  (1970)  gives  an  historical  re- 
view). Guldberg  introduced  this  distribution  in  1934,  Lundberg 
first  used  it  in  1940  as  a model  for  accident  proneness  and 
Arbous  and  Kerrich  (1951)  expanded  G uldberg' s work  and  fitted 
the  distribution  to  bivariate  data  related  to  accidents  in  in- 
dustrial settings.  In  addition  to  their  theoretical  contribu- 
tions Bates  and  Neyman  fitted  the  distribution  to  numerous  bi- 
variate data  sets  related  to  diseases  in  industrial  workers. 

In  1954  Arbous  and  Sichel  fitted  it  to  shift-workeT  absenteeism 
data  for  adjacent  time  periods.  Youngs,  Geisler,  and  Brown  in 
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1955  studied  the  conditional  distribution  of  this  bivariate 
r.v.  and  showed  how  it  could  be  used  for  the  prediction  of 
demand  for  aircraft  spare  parts.  In  1961  Edwards  andGurland 
generalized  (2.6)  and  compared  the  fits  obtained  from  the  two 
distributions.  The  regression  function  for  their  distribution 
is  linear  also.  Subrahmaniam  (1966)  and  Subrahmaniam  and 
Subrahmaniam  (1973)  also  studied  this  latter  bnb  distribution. 

Rjr  purposes  of  comparison  we  choose  to  call  the  bnb  dis- 
tribution in  (2.6)  the  G uldberg-Bates -Neyman  model  with  para- 
meters a,  6 and  v and  to  designate  it  G -B-N(ot,9,v)  . 

Certain  data  sets  do  not  exhibit  empirical  regressions 
which  are  linear  nor  do  some  data  sets  show  positive  correla- 
tion and  so  it  is  natural,  for  these  cases,  to  work  with  a 
bivariate  probability  function  which  allows  for  nonlinear  re- 
gressions or  negative  correlations  or  both.  The  classical 
Bates  and  Neyman  paper  exhibited  empirical  data  best  fitted 
by  regression  curves  which  were  obviously  nonlinear,  and,  con- 
sequently, their  results  were  not  entirely  satisfactory.  fur- 
thermore, we  show  a new  bivariate  data  set  related  to  aircraft 
flight  aborts  which  has  a negative  sample  correlation  coefficient. 

In  the  next  section  we  discuss  a new  bnb  distribution  which 
admits  of  nonlinear  regressions  and  negative  correlations  and 
derive  several  of  its  properties. 
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2 . 2 A New  Bivariate  Negative  Binomial  (and  Gamma)  Distribution, 
l£a  Convolution,  and  Some  Properties 

The  purpose  of  this  section  is  to  show  how  a new  bnb  dis- 
tribution may  be  obtained,  by  the  process  of  convolution,  from 
a certain  bivariate  geometric  distribution.  A number  of  proper- 
ties such  as  the  moments  to  order  two,  the  regression  function, 
a recursive  formula  for  the  cell  probabilities,  and  the  likeli- 
hood equations  are  obtained  for  certain  special  cases. 

Paulson  and  Ifcpuluri  (1972)  showed  that  the  bivariate  r.v. 
(X  ,Y)  , where  each  element  in  the  pair  is  defined  on  the  non- 
negative integers,  has  a bivariate  geometric  distribution  if 
its  characteristic  function,  ♦ satisfies  the  character- 

istic-functional equation 


♦ CD  = <KD  E[  ♦(TV)]  (2.13) 

where 

T “ ^tl»t2^»  ♦CT)  " ♦iC1!*0)  ♦2^0,t2^» 

♦iCtj.O)  = [1+jSp  (l-e1*!)]’1, 

♦ 2(0.t2)  " t1+  bq  U-*1*2)]’1. 

and  V is  a 2x2  matrix- valued  r.v.  having  values  in  the  set 

{ [o  o]  * [o  o]  » [o  l]  * [o  l]  } with  Probabilities  a,b,c 
and  d,  respectively.  Here  0 <p  <1,  0 <q  <1,  a+b+c+d-1, 


b+d  <1  and  c+d<l.  Thus  the  characteristic  functional  equation 
can  be  rewritten  as 


♦ (T)  - H>1(t1,0)  H»2(0,t2)[a+b<Ktlt0)  + c$(0,t2)  + dp  (T)  ] (2.14) 
and  it  is  easy  to  show  that 


4>  Ctx  ,0) 
<t>(0,t2) 


(a+c)^1 

1-  fb  + dT'i'Y 

(a+b)4»2 
r-  (c+d)  p2 


[l+01(l-eitl)]”1, 

[l+02(l-eit2)]-1 


(2. IS) 


whe  re 


6X  * P/[(1-P) (a+c)] , 02  = q/[(l-q) (a+b)]  (2.16) 

and  the  arguments  of  i^Ct.^,0)  and  \p2(0,t2)  have  been  suppressed 
(and  will  be  in  the  sequel).  Comparing  (2.15)  and  (2.2)  we  see 
that  the  marginals  are  geometric. 

In  (2.13)  the  characteristic  function  ip  (T)  = i|^(t^,0)- 
i^2(0,t2)  corresponds  to  a bivariate  geometric  distribution  with 
independent  geometric  marginals  and  as  Block  (1975)  points  out, 
(2.13)  gives  the  characteristic- functional  equation  of  the  bi- 
variate random  variable 

N1  N2 

(X,Y)  - ( Z I Y.).  (2.17) 

i=l  1 i*l  1 

The  pair  (N^,N2)  is  a certain  bivariate  geometric  distribution 
and  is  independent  of  (X^,Y^),  i = l,2,3,...,  which  are  indepen- 
dent and  identically  distributed  (i.i.d.)  random  variables  with 
characteristic  function  ip (T) . Thus  (2.13)  corresponds  to  a 
special  kind  of  bivariate  geometric  compounding  of  the  distri- 
bution with  characteristic  function  ip (T) . It  is  also  clear  that 
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if  <Ji(T)  has  geometric  marginals,  then  <J>(T)  will  have  geometric 
marginals.  This  is  even  more  apparent  from  (2.17)  since  uni- 
variate random  geometric  sums  of  i.i.d.  geometric  random  vari- 
ables are  geometric. 

Paulson  (1973)  has  shown  that  (2.13)  also  characterizes  a 
bivariate  exponential  distribution  where  iHT)  is  the  character- 
istic function  of  independent  exponential  r.v.'s.  In  addition 
he  forms,  in  a way  to  be  described  here  for  a bnb  distribution, 
a bivariate  gamma  distribution.  Certain  properties  to  be  de- 
rived for  this  bnb  distribution  will  be  extended  to  his  bivari- 
ate gamma  distribution. 

Paulson  and  Ifcpuluri  obtained  the  moments  of  the  distribu- 
tion of  (X  ,Y)  in  (2.13)  to  order  two  and  showed  that  the  corre- 
lation coefficient  has  values  in  the  interval  -0.25  _<  p <1. 
They  also  presented  recursive  formulae  for  determining  the 
probability  function. 

Clark  (1972)  obtained  a closed  form  representation  for 
the  bivariate  geometric  distribution  characterized  by  (2.13) 
for  the  special  case  b=c=0  and  extended  it  to  a bnb  distribu- 
tion. Next  we  summarize  that  development.  (By  defining  a 
multivariate  analogue  of  the  characteristic- functional  equation 
(2.13),  Clark  also  constructed  a multivariate  geometric  distri- 
bution and  extended  it  to  a multivariate  negative  binomial  dis- 
tribution; Appendix  A shows  the  unpublished  derivation.) 


""  r' 
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We  consider  the  case  b=c=0  in  (2.14)  and  so 

♦ (T)  = 4'1'l'2  (a+d<ji  (T) ) . (2.18) 

Solving  for  d> (T)  leads  to 

4>(T)  = a^^l-dij,^]"1  (2.19) 

and  upon  expansion 

♦ (T)  = a^»1^2  [l+di^1^2  + (d\p1ip2 ) 2 + ...]  . (2.20) 

The  inverse  transform  of  <j>(T)  , that  is,  the  probability  func- 
tion, say,  g1(x,y),  may  be  obtained  termwise  from  (2.20)  since 
the  resultant  series  converges  uniformly  and  absolutely  for  all 
x,y-0 , 1 , 2 , . . . , (Titchmarsh(1964));  we  have 

oo  t m 

g (X,y)  = a(l-p)px(l-q)qy  l (X+xj)(y^)[d(l-p)(l-q)]J  (2.21) 

1 j=o  x y 

where  x,y=0 , 1 ,2 , . . . . Expanding  the  combinatorial  terms  gives 

gj/x.y)  a a(l-p)px(l-q)qy  F(x+l,y+l;  1;  d(l-p)(l-q))  (2.22) 

where  F(a,b;c;z)  is  the  Gaussian  hypergeometric  series  given 

by 

- Ca) . Cb)  . j 

F(a,b;c;z)  « 1 + £ — ■ jr  (2.23) 

j = l 

and  the  term  (n)^  is  defined  by 

(n)j  * = "(n+l)  (n+2)  . . . (n+j-1)  . (2.24) 
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Taking  the  v-fold  convolution  of  gjO.y)  with  itself 
yields  a bnb  distribution  which  is  denoted  by  gv(x,y).  The 
characteristic  function  of  gv(x,y)  is 


4>V(T)  = [4>(T)]V  = (a^1i/»2)v  [l-d^1^2] 
= (a4»1ii>2)v  [l+v(d^1^2)  + I (di^^)2 


and  in  the  same  manner  as  before  we  obtain 


-v 


...]  (2.25) 


gv(x,y)  = avh1(x)h2(y)  F(x+v ,y+v ;v ; z) 


(2.26) 


where 


hl(x)  = (^Hl-plV. 


h2(y)  = (71)(l.q)V, 


(2.27) 


z = d(l-p) (1-q) 

and  x,y=0 , 1 , 2 , . . . . Of  course,  v=l  leads  to  (2.22).  It  can 
be  shown  that 


♦v(tlf0)  » [1+61 (l-eitl) ] "V 


and 


ov(o.t2)  = [i*e,u-eit2)]'v 


(2.28) 


! 


where  9^,  and  02  are  defined  as  in  (2.16)  and  so  the  marginals 
of  (X,Y)  are  negative  binomial. 





Clark  obtained  the  moments  to  order  two  of  this  bnb  dis- 
tribution; the  remainder  of  his  work  was  limited  to  numerical 
investigations  of  the  bivariate  geometric  distribution  charac- 
terized by  (2.14)  with  b^O  and  c/0 , and  he  showed  figures  of 
the  probability  surface  and  the  regression  function  as  it  de- 
pends upon  p,  the  correlation  coefficient.  The  probability 
surface  was  computed  by  using  the  recursive  formulae  given  by 
Paulson  and  l^puluri  and  then  the  regression  function  was  com- 
puted by  using  the  definition  of  a conditional  mean,  that  is, 
E[Y|x]  = £ y Pr(y|x),  where  the  summation  is  taken  over  all 
non-negative  y.  His  numerical  results  showed  that  the  proba- 
bility function  admits  of  nonlinear  regression  functions  with 
either  positive  or  negative  correlations. 

Next  we  generalize  Clark's  work  for  the  cases  b and  c not 
equal  zero  and  show,  among  other  things,  an  analytical  deriva- 
tion of  the  regression  function  for  some  special  cases. 

We  construct  in  a way  parallel  to  Paulson's  (1973)  deri- 
vation for  a bivariate  gamma  distribution  a bnb  distribution. 
From  (2.14) 

$(T)  * ^1^2  (a  + + c<f>2  + d <#>  (T)  ] (2.29) 

or 

♦(T)  ’ (l-d»^2)ta  * b*l  * c*2>  C2.30) 

and  convolving  as  in  (2.25)  yields 
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♦vm  ■ * I>*1  * «2]v.  (2.31) 

From  (2.29)  on  we  write  and  <J>2  for  <p(t^,0)  and  4>(0,t2), 
respectively.  Recalling  from  (2.13)  that  ^ and  i P2  are  func- 
tions of  p and  q respectively,  we  choose  to  designate  the  par- 
ticular bnb  distribution  which  results  here  as  the  BNB(a,b ,c,p 
distribution  and  to  label  the  probability  function  as  fv(x,y). 
It  is  relativel)  easy  to  show  from  (2.31)  that  the  marginals 
are  negative  binomial.  The  following  theorems  and  results  are 
now  presented. 


Theorem  1:  The  inverse  transform  of  4>V(T)  defined  by  (2. 

with  a,b,c,d  as  probabilities  and  a+b+c+d=l,  b + d<l,  c+d<l  and 
integer  is 


fv(x,y) 


l 

a.e.Y 


v!  _a-v1_3_yA 

oTSTyT  3 b C *1.0 


* gv(x,y)  * 


r «i 

2 -Y 


(2.32) 


where  £ runs  over  all  cx,B,y  > 0 such  that  a ♦ B + y * v , 
x,y=0 , 1 , 2 , . . . , gv(x,y)  is  the  probability  function  for  the 
BNB(a,0,0,p,q,v)  of  equation  (2.26)  and 


>1.3  * (m-)8CXTl)(T^7)X' 


*2.y  - (^YfyT1)(I^T)y 


(2.33) 


x 

The  operator  * for  convolution  over  x is  defined  for  two  func- 
tions h^x.y)  and  h2(x,y)  by 

X X 

h.  * h2  - l h. (£,y)h2(x-£,y) . 

1 c C-0  1 L 


31) 

v 


(2.34) 


I 

[ 


16 


The  operator  * is  defined  similarly. 


Proof:  Prom  (2.31)  and  using  the  multinomial  expansion 


gives 


*1*2 


*v<T>  “ ^l-diKl'J  £ a’.g.'y'.  a tb  + i)S(c^2)Y  C2-3S) 

i i a,t > ,y 

where  the  £ is  over  all  a,8,y  _>  0 such  that  a ♦ 6 + y = v . 
Recall  from  (2.15)  that  4^  = [l  + flj  (l-eit:l)  ] ” 1 and 
<t>2  * [l-^02  C1_elt2)]  1 and  changing  to  z-transforms  by  letting 
e**!  * r 1 and  e^t^  = s * we  can  write 


t+§y 


and  <p2s 


^7  e2 

S'TTW2 


where  <J>lr  and  d> 2 s are  the  corresponding  transformations.  Taking 
the  inverse  z-transform,  ^_1,  of  <J>jr  and  <J>^s  we  get  (from  the 
z-transform  pair  number  24,  Jury  (1964)) 

f lC*!r>  ■ *1.6  - ♦l.o'1’ 


f'M.’  ■ *2.y  * (T^)Y(yT1)(T^I)y-  Y^- 

Except  for  the  constant  aV,  the  inverse  transform  of  (^  ^ )v 
is  gv(x,y)  in  equation  (2.26)  so  (2.35)  may  readily  be  inverted 
(after  changing  to  z-transforms  throughout)  by  first  holding  s, 
say,  constant  and  inverting  with  respect  to  r and  then  completing 


t 
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the  inversion  by  inverting  with  respect  to  s.  Thus  we  arrive 
at  (2.32).  This  equation  is  analogous  to  Paulson's  (1973)  re- 
sult for  a bivariate  gamma  distribution  (aa  in  equation  (26) 
of  that  paper  should  be  replaced  with  aa  v) . Directly  from 
Theorem  1 we  have 

Corollary  1:  For  v=l  and  a^O , the  probability  function 

for  the  BNB(a,b,c,p,q,l)  distribution  is 

f^x.y)  ■ gjfx.y)  * |(T^')crr5Y)x  * M*./) 

* f(T47>^)y  1 M*-/)  (2.36) 

where  g1(x,y)  is  the  probability  function  for  the  BNB  (a  ,0  ,0  ,p  ,q  , 1) 
of  (2.21).  This  is  a closed  form  representation  of  Paulson  and 
Uppuluri's  bivariate  geometric  distribution. 

Obviously,  the  utility  of  fv(x,y)  in  (2.32)  is  limited  by 
the  integral  requirement  for  v and  so  we  seek  a representation 
for  v real  valued.  Except  for  the  case  c=0  our  attempts  to  de- 
rive a general  representation  have  been  unsuccessful.  Next  we 
give  the  results  for  c=0  (or  for  b=0  by  symmetry).  Also,  we 
show  how  these  results  are  directly  applicable  to  Paulson's  bi- 
variate gamma  distribution. 

Theorem  2:  Rjr  the  BNB (a ,b ,0 ,p ,q ,v)  distribution  with 

v>0 , 

X oo 

fv(x,y)  - gvCx,y)  ♦ gv(x,y)  * J h(b ,k) (l-p)k(x+^'1)px  (2.37) 


r 

f-  f- 


I 
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where  gv(x,y)  is  the  BNB (a  ,0 , 0 ,p  ,q , v)  of  f 2 . 26 ) and 

h(b,k)  = [-(b+d)]k  X OOC)n-  (2.38) 

n=0 

Proof:  In  (2.31)  we  write  0 (c=0)  for  the  characteristic 

function  when  c=0  and  by  using  (2.15)  it  follows  that 


a*l*2  ,», 


b*. 


Vc=0)  ~ ci^ip1ip2)  [1+rnwnq-i 


(2.39) 


The  leading  term  in  this  expression  is  the  characteristic  func- 
tion for  the  case  b = c=0  and  for  it  we  write  <f>v(b=c=0).  Re- 
writing the  second  term  and  expanding  in  an  infinite  series 
gives 

$v(c  = 0)  = 4>v(b=c=0){r-nr^T^)v 

(JO 

- ^ (b=c*0)  (1+  I [ l (J.jJ  («)  (b+d)k"ndnJ  (-^1)k> 
v k=l  n=0  K n n 1 


= <p  (b=c*0){l+  l h(b  ,k)  i|>,k } 
v k=l  1 


(2.40) 


where  h(b,k)  is  defined  in  (2.38).  Since  the  infinite  series 
in  (2.40)  is  uniformly  convergent  we  can  invert  termwise  to  get 
(2.37).  It  is  easily  verified  that  b=0  leads  to  gv(x,y)  in 
(2.26)  and  v=l  gives  f1(x,y)  in  Corollary  1 with  c=0. 

As  stated  before,  equation  (2.13)  characterizes  a bivariate 
exponential  distribution  and  for  that  case  (tj  ,0)  (°  “ 

[l-91it1] '1[l-e2it2]"1  where  0j  and  02  are  different  from  those 
in  (2.16).  From  the  characterization,  Paulson  (1973)  formed  a 
bivariate  gamma  distribution  and  for  the  special  case  of  c=0  a 
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corresponding  result  to  Theorem  2 is 

x oo  k-1  -x/9 

fv(x,y)  * gv(x,y)  + gv(x,y)  * ^ h(b,k)(^-)k  - — — - (2.41) 


gv(x,y) 


(Ie^T^2),i(V'l,''X/6l-y/92  (*•«) 


h(b,k)  is  the  same  as  in  (2.38),  and  Iv_^(2(q-*-^  )**)  is  the  mod- 
ified Bessel  function  of  the  first  kind  and  order  v-1.  That 


this  result  is  true  follows  directly  from  (2.40)  since  ^ in 
this  situation  is  [l-0^it^]-1.  lor  v=l  and  c=0,  (2.41)  checks 
with  a result  by  Kohberger  (1975).  In  a slightly  different 
form  (2.42)  is  the  bivariate  gamma  distribution  obtained  im- 
plicitly by  Wicksell  (1933)  and  explicitly  by  Kibble  (1941). 
The  bivariate  exponential  distribution  defined  by  (2.13)  will 
be  discussed  in  more  detail  in  Part  4.  We  continue 'with  some 
more  properties  of  the  BNB(a ,b ,c,p ,q,v)  distribution. 

By  using  the  characteristic  function  4>v (T)  in  (2.31)  and 
the  usual  differentiation  techniques  there  follows,  after  sev- 
eral tedious  operations. 


Theorem  3;  The  mean  vector  and  covariance  matrix  for  the 
BNB(a,b ,c,p ,q,v)  distribution  are 


\ 

I 


r- 

i 


I 
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r-' 


and 


vejCi-0!)  V(^-y.8l02 

''(ajld':Jei92  '’e2ci-e2) 


(2.43b) 


We  digress  briefly  at  this  point  and  give  a method  of  moments 

parameter  estimation  technique  for  the  special  case  b=c=0. 

1.  In  an<*  °y  set  91  = x/v  and  02=y/v  and  take  the  pro- 
2 2 

duct  of  o£  and  Oy  to  be  equal  to  the  product  of  the  sample 
variances.  There  results  a quadratic  function  in  v: 


-V 

(1  - — — -^)  v2  ♦ (x+y)v  ♦ xy  * 0 
xy 


(2.43c) 


.2  2 


and  for  s^Sy  > xy,  which  is  expected  if  the  marginals  are 
approximately  negative  binomial,  there  is  exactly  one  posi- 
tive root  that  we  can  take  as  our  estimate  of  v. 

2.  In  c^y  substitute  0^*x/v  and  02=y/v  and  set  a^y  eclual 
to  the  sample  covariance,  s^y.  Solving  for  the  parameter  a, 
we  have 

, VSXY 


(2.43d) 


xy 

3.  from  (2.16),  9-j^p/  [(l-p)a]  and  so 

a0. 


...  1 
I+a’e. 


ax/v 

1+ax/v 


(2 . 43e) 


4.  In  a like  manner. 


1+ay/v 


(2 . 43f) 


I 


V.. 


■ *'  sty 
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The  next  theorem  establishes  that  the  regression  function 
for  the  bivariate  geometric  distribution  in  Corollary  1 is  non- 
linear. 


Theorem  4:  For  the  BNB  (a  ,b  ,c  ,p  ,q  ,1)  distribution  (bivar- 

iate geometric  of  Paulson  and  Uppuluri)  where  b^O  the  regression 
function  is 

Enr|x]  - A It1*  <s-  F)kX*1'  (2 .44) 

where 


m * p + (1-p) (a+c)  , A=cm/ [(a+c) (a+b) ] , 
k = m/ [m  + b(l-p)]  and  d = 1-a-b-c.  (2.45) 


Proof:  For  the  conditional  mean  E[Y|x],  which  is  a 

function  of  x,  take  as  a generating  function  the  z-transform 


^(E[Y|x])  = g(z)  = l z"X  E [ Y | x ] . 


(2.46) 


By  using  the  definition  of  E[Y|x]  we  can  write 

» * U+M  x 

g(z)  = (1  + 0,)  l l [ -=g  ]X  yfjCy.x)  (2.47) 
1 x=0  y=0  Z01  1 

1 ®i  x 

since  the  marginal  density  of  X is  Pr[X=x]  * (jyg— ) (^-g-  -)  . 
This  last  equation  can  be  written  as 

1+6, 


gU)  - (1‘ep  & Pl-jji.xilx  . J 


(2.48) 
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where  P[-,’]  is  the  probability  generating  function  of  the 
bivariate  pair  (X,Y) . The  probability  generating  function 
can  be  obtained  from  4>V(T)  by  substituting  r=eltl  and  s=ext2 
and  after  some  computation  and  simplification  we  get 

g(z)  = (m+b  cr-p) ) + rz-'ki  Cz-rr}  (2<49) 

where  m,  k and  A are  defined  in  (2.45).  Inverting  g(z)  (Jury 
(1964))  gives  (2.44).  This  type  of  a nonlinear  regression 
function  is  sometimes  called  an  exponential  regression  function. 
In  a direct  way  using  the  definition  of  E[Y|x]  we  obtain 

Theorem  5:  For  the  BNB (a ,0 ,0 ,p ,q ,v)  distribution 

EtYM  = (2.50) 

where  d=l-a. 

The  following  theorem  is  very  useful  for  computations. 

Theorem  6:  For  the  BNB  (a  ,0 ,0  ,p  ,q , v)  distribution  the 

probability  function  can  be  computed  with  the  recursive  formula 

gv(x,y+l)  - (yVl^l.-zyfCx^vjg^x^j-p^x-Dg^x-l.y)],  (2.51) 
x > 1,  y > 0 and  z - d(l-p)(l-q).  Coupled  with 

gv(0,0)  - 


the  probability  function  is  determined. 
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Proof:  By  using  the  following  relation  for  F(a,b;c;z) 

and  two  of  its  contiguous  functions  (15.2.18,  Abramowitz  and 
Stegun  (1964) ) : 

(c-a-b) F (a  ,b ;c  ; z)  - (c-a) F (a- 1 ,b ;c ; z) 

♦ b(l-z)F(a,b+l;c;z)  = 0, 

(2.26)  can  be  written  as 

gv(x,y)  = avh1(x)h2(y)  (^~r[xF(v+x-l,v-»y;v;z) 

+ (v+y) (1-z) F(v+x,v+y+l  ;v ; z) ] } 


and  (2.51)  follows  easily  with  some  elementary  operations. 

The  cell  probability  at  (0,0)  is  obvious  from  (2.26)  and  (2.23). 
An  important  result  for  parameter  estimation  is  next. 


Theorem  7:  For  the  BNB  (a  ,0  ,0  ,p  ,q  ,v)  distribution  and  if 

v is  known  the  likelihood  equations  for  a random  sample  of  size 
n are 


A±°& 

3 a 


L. 


(2.52a) 


•blog  L. 
3 P ' 


- K 


0 


(2.52b) 


±l2&-k:  ikalX  ♦ y - R , o 

3<1  L * 


(2.52c) 


where  L is  the  likelihood  function, 
for  the  marginal  distributions 


x and  y are  the  sample 

n-1  t „ r'r*1i8',0c,y*1) 

’ » X.y  Vc  q ’8v(x,y)  • 


means 
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n is  the  number  of  observations  for  which  X=x,  Y=y,  and 
xy 

gv(x,y)  is  the  probability  function  in  (2.26). 

Proof:  If  the  probability  function  in  (2.26)  is  differ- 

entiated with  respect  to  the  parameters  a,  p and  q the  follow- 
ing differential-difference  equations  result: 

“X5 + 5)gv(x»y)  - Z(Iqi)gv(X’y+1)  (2.53a) 

*g,.(x,y)  vv  1 v+1 

~ VP = CP  + T^P)gv(x»y)  • Fp(V)gv(x’y+1)  (2.53b) 

"b  gv(x,y) 


• p'  °v  v ' 1-p v q 

(q  + T^q^v(x'>r)  " ITq(ZVL^v(x'y+1)  * (2.5 


(2.53c) 


These  equations  follow  by  using  (15.2.1,  Abramowitz  and  Stegun) 

>FO,b;c;i;)  - (a+l,b+l;c+l;z) 

^ z c 

and  (exercise  1,  page  296,  Whittaker  and  Watson  (1965)) 


F(a,b+l;c;z)  - F(a,b;c;z)  * |^-F(a+1  ,b  + l ;c+l ; z)  . 

The  log  likelihood  function,  log  L,  for  a random  sample 
of  size  n is  l nxy  log  gXJ(x,y)  and  so 


x,y 


blog  L y n 1 *gU(X>y) 

ia  x^y  xy  gv (x,y)  * a 


Using  (2.53a)  and  a few  simple  operations  leads  to  (2.52a). 
Similarly,  (2.52b)  and  (2.52c)  obtain.  We  were  unable  to  get 

>gv(x,y) 

a differential-difference  equation  involving  — — , v as- 


b v 


sumed  unknown. 
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From  (2.52)  it  is  clear  that 


P 


P.:Rly  = V 

q ’ a 


(2.54) 


and  these  relationships  are  very  useful  in  estimating  the  para 
meters  via  the  method  of  maximum  likelihood.  For  v known  or 
not,  the  conditions  in  (2.52)  are  necessary  for  a maximum  like 
lihood  solution  to  the  likelihood  function  for  b=c=0.  There- 
fore, (2.54)  can  be  used  to  reduce  the  dimensionality  of  the 
unknown  parameter  space  from  four  if  v is  unknown,  to  two  by 
taking,  say,  p = qx/ [qx+ (l-q)y ] and  a = vq/[(l-q)y].  We  have 
used  a nonlinear  optimization  computer  program  (Cross  (1970)) 
to  solve  for  the  parameter  estimates  and  the  dimensionality 
reduction  permits  extremely  shorter  running  times. 

Next  we  show  a corresponding  result  to  Theorem  7 for 
Paulson’s  bivariate  gamma  distribution.  For  the  distribution 
defined  in  (2.42)  and  if  v is  known  the  likelihood  equations 
for  a random  sample  of  size  n are: 


>L°g-L. 

f 

V 

- S = 

1 a * 

a 

blog  L.  * 
* B1  ’ 

X 

91 

V+l  c- 

~1  s " 

"blog  L. 

J. 

- - 5 - 

t>q7  ' 

. 07 

2 3 

0 

0 


(2.55) 


where  L,  x and  y are  as  in  Theorem  7,  a=l-d,  \ I*/I, 

x,\ 

I*-(z/2) ^ and  z-2[(l-a)xy/(8j62)p-  The 


result  follows,  after  some  lengthy  but  straightforward  computa 
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tions,  by  taking  L = n g (x,y)  where  g (x,y)  is  in  (2.42) 

x,y 

and  forming 

>1o«  L y 1 

a*  x,y 

for  Xe{a, 0 p 9 2 ^ ^ is  obvious  that  necessary  conditions  for 
the  parameters  are 


We  show  in  a later  section  the  bnb  of  this  section  fitted 
to  some  data.  First  we  introduce  another  bnb  distribution 
which  has  certain  desirable  properties. 

2 . 3 A New  Bivariate  Negative  Binomial  Distribution,  Via  a 

Trivariate  Reduction,  and  Some  Properties 

The  two  previously  discussed  bivariate  negative  binomial 
distributions  have  marginals  which  are  negative  binomial  with 
parameters  0^,  i-1,2,  and  common  parameter  v.  In  this  section 
we  introduce  another  bivariate  negative  binomial  distribution 
whose  marginals  have  parameters  , i=l,2,  and  common  parameter 
6.  Data  are  shown  later  for  which  this  latter  model  seems  more 
appropriate. 

We  construct  via  reduction  of  a certain  trivariate  nega- 
tive binomial  distribution  with  independent  marginals  a bivari- 
ate negative  binomial  distribution.  Mardia  (1970)  refers  to 
this  as  a trivariate  reduction;  Holgate  (1964)  used  this  tech- 
nique to  construct  a bivariate  Poisson  distribution  and  Arnold 
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(1967)  generalized  the  procedure. 

Theorem  8:  Let  X^,  X2  , X^  be  independent  negative  bi 

nomial  r.v.'s  with  common  parameter  0 and  individual  para- 
meters Vj,  and  v^,  respectively.  Then  the  probability 

function  of  X = Xj  + X2,  Y = X2  + Xj  is  given  by 


where 


. , ...  1 . Wv3,  0 ,y  ? „ 

h(x,y)  - (1  + 0)  ^1  + 6-^  ^ ^ w 

W K 


v2+x-w-l  v3+y-x+w-l 
r L ~)(  1 ( .JL  -) w 

^ V-I.f  ' ^ V - v4- ur  7 vi+0' 


x-w 


y-x+w 


(2.57) 


k = 


0 , if  x < y 
x-y,  if  x > y. 


Proof:  The  joint  distribution  of  , X2 


X3  is 


x- 


,,  . . i . vi+v2+v3  * 

fCx1,x2,xJ)  - C1+e)  ^ ^1+6^ 

and  by  taking  the  transformation  of  variables  X = Xj  + X2 , 

Y * X2  ♦ Xj  and  W = X^  it  follows  in  the  usual  way  that  the 
joint  distribution  of  (X,Y)  is  (2.57). 

This  bivariate  negative  binomial  distribution  is  designated 
BNB-TR(0 ,Vj ,v2 ,v3)  and  from  the  defining  relations  its  marginals 
are  negative  binomial;  X has  parameters  9 and  + v2  and  Y has 
parameters  0 and  v2  ♦ (Johnson  and  Kotz).  The  marginal  means 
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and  variances  are  thus  known  and  the  covariance  of  (X,Y)  is 

cov(X,Y)  = cov(X1  + X2,  X2  + X3)  = c|  (2.58) 

Therefore,  p = v2/[(v2+v2)  ^v2+v3^  ] 2 an<^  we  see  that  p is  re- 
stricted to  nonnegative  values. 


Theorem  9:  The  regression  function  for  the  bivariate  r.v. 

in  (2.5  7)  is 


E[Y|xl  = v3e  ♦ .v-£_.  (2.59) 

JL  L 

Proof:  Given  X=x,  the  r.v.  Y |x  has  expectation  E[X2|x] 

+ v since  X^  and  X are  independent  and  so  we  desire  the  dis- 
tribution of  X2|x.  By  writing  the  joint  distribution  of  Xj 
and  X2  and  transforming  to  new  variables  by  letting  X = Xj  + X2 , 
*2  “ *2*  we  °htain  for  the  joint  distribution  of  (X,X2) , 


, ,Vx-x2-\f',2*x2-1w  1 ,VV2,  6 ,x 

f(x,x2)  - ( x.Xz  )(  Hiry)  (iTj)  • 


0 <_  x2  £ x. 


(2.60) 


From  definitions  it  follows  that 


f(x,|x) 


u.+x-x9-l  v7+x,-l 

( 1 2 ) ( 2 2 ) 

v x-x2  x2 

V. +V-+X-1 

( V ) 


, 0<x2£x 


(2.61) 


and 


x 


I 
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pare  with  the  Guldberg-Bates-Neyman  model  in  fitting  bivariate 
data  related  to  shift  worker  absenteeism  and  disease  data  for 
industrial  workers. 


2 . 4 Bivariate  Negative  Binomial  Distributions  Fitted  to  Data 
In  this  section  we  fit  two  data  sets  from  the  literature 
with  the  three  previously  discussed  bnb  distributions.  Our 
objective  in  using  these  data  is  to  illustrate  certain  aspects 
of  the  distributions.  The  first  data  set  is  given  by  Arbous 
and  Sichel  (1954)  and  concerns  absenteeism  for  248  shift  workers 
in  two  adjacent  yearly  time  periods  and  the  second  one  is  due 
to  Bates  and  Neyman  (1952)  and  shows  the  number  of  cases  of  in- 
capacity suffered,  per  individual,  during  a common  time  period 
and  due  to  two  diseases;  the  sample  size  is  1286.  For  the  ab- 
senteeism data  we  show  that  either  the  G-B-N (a , 9 ,v)  model  in 
(2.6)  or  the  BNB (a , 0 , 0 ,p ,q , v)  model  in  (2.26)  fit  the  bivariate 
data  reasonably  well  but  that  the  regression  function  from  the 
latter  model  describes  better  the  observed  conditional  means. 
None  of  the  three  models  adequately  describe  the  joint  distri- 
bution of  the  disease  data  even  though  the  marginals  are  accep- 
tably fitted  by  a univariate  negative  binomial  distribution.  It 
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will  be  shown,  howeve r,  that  the  two  new  bnb  distributions  do 
describe  certain  properties  of  the  disease  data  better  than 
the  G- B-N (a , 9 , v)  distribution,  though. 

Parameter  estimates  for  the  various  models  are  by  the 
method  of  maximum  likelihood  estimation  (MLE) . Bates  and 
Neyman  show  how  to  solve  the  likelihood  equations  for  the 
G-B-N(ct,9,v)  distribution;  for  the  BNB  (a  ,b  , c ,p  ,q  ,v)  and  BNB-TR 
(0,Vj,V2,Vj)  distributions  we  use  a computer  optimization  pro- 
gram (Cross)  directly  on  the  log  likelihood  function  to  obtain 
our  estimates.  For  the  parameters  b and  c defined  to  be  zero 
in  the  BNB  (a  ,b  ,c  ,p  ,q  ,v)  model  we  use  the  necessary  conditions 
in  (2.54)  to  reduce  the  dimensionality  of  the  parameter  space. 

Corresponding  to  the  observed  pairs  (x^y^),  i=l,2,...,n, 
representing  a random  sample  from  the  unknown  probability  func- 
tion f(x,y),  we  wish  to  test  the  hypothesis  HQ  : f (x  ,y)  =*fQ  (x  ,y)  , 

where  fg(x,y)  is  specified  to  be  one  of  the  referenced  bnb  dis- 
tributions.  The  x test  is  used  as  a goodness-of-fit  test  of 
Hq;  we  point  out  that  the  problems  with  grouping  cells  which 
are  generally  associated  with  this  test  in  univariate  settings 
are  even  more  dramatic  for  bivariate  cases.  To  illustrate, 

Table  10  shows  two  independent  groupings  for  a data  set  to  be 
described  in  the  next  part.  Before  commenting  on  the  table, 
certain  remarks  are  required. 

Following  the  practice  of  Bates  and  Neyman  in  their  paper, 
cells  which  have  expected  frequencies  less  than  three  are 
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grouped.  In  Table  10  and  similar  ones  to  follow,  the  dashed 

lines  indicate  the  boundary  of  the  particular  cells;  heavy 

lines  indicate  the  grouping  adopted  for  the  application  of  the 
2 

X test.  Three  numbers  are  shown  in  each  cell:  the  observed 

frequencies  are  shown  in  the  upper  left  corner  of  particular 

cells  and  the  decimal  numbers  are  the  expected  frequencies  on 

2 

the  left  and  the  contributions  to  x on  the  right.  If  several 

2 

adjoining  cells  are  grouped  then  the  expected  frequency  and  x 

values  shown  are  for  the  entire  group.  The  P value  given  is 

2 

the  probability  of  obtaining  a value  of  x exceeding  the  com- 
puted amount  assuming  is  true. 

Table  10  clearly  shows  how  the  probability  P is  affected 
by  different  groupings . Unlike  the  univariate  situation  we 
have  two  directions  to  contend  with  here  for  grouping  and  it 
is  not  obvious  how  to  proceed.  This  is  pointed  out  to  empha- 
size the  need  for  an  alternate  goodness -of - fit  test  for  bivari- 
ate data  (and  in  general  multivariate  data) , one  perhaps  being 

independent  of  any  grouping.  For  lack  of  a better  test  we  re- 
2 

sort  to  the  x test.  In  every  fit  to  be  shown  herein  the 
groupings  are  completely  independent  of  the  observations  and  the 
findings  are  from  a single  attempt  at  grouping  the  expected  fre- 
quencies for  the  minimum  value  of  three. 

Table  1 shows  certain  summary  results  for  all  of  the  data 
to  be  analyzed.  The  first  column  identifies  the  data,  column 
2 gives  the  sample  size  and  correlation,  column  3 specifies  the 
marginal  random  variables  and  shows  the  associated  sample  means 
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and  variances,  columns  4 and  5 show  the  parameters  of  the  uni- 
variate negative  binomial  fitted  to  the  marginals  and  the  assoc- 
2 

iated  x values,  degrees  of  freedom  and  probability  levels,  re- 
spectively. Results  from  this  table  will  be  presented  along 
with  a discussion  of  the  individual  data  sets.  The  absenteeism 
data  is  examined  first. 

One  of  the  objectives  of  the  Arbous-Sichel  paper  was  to 
extend  the  notion  of  accident-proneness  introduced  by  earlier 
workers  (see  Kemp  (1970)  for  a history)  to  absence -proneness 
in  shift  workers.  Table  2 shows  the  observed  and  expected  cell 
frequencies  for  the  number  of  absences  in  two  adjacent  yearly 
time  periods  (1947  and  1948)  for  248  workers;  the  expectations 
are  from  the  G-B-N(l,0,v)  symmetric  model  (ot=l)  . 

From  Table  1 we  see  that  the  marginal  distributions  are 
fitted  rather  nicely  by  the  univariate  negative  binomial  and 
coupled  with  the  fairly  large  sample  correlation  coefficient, 
it  seems  reasonable  to  expect  a bnb  distribution  to  adequately 
describe  the  data. 

Table  3 shows  the  expected  cell  frequencies  from  the  BNB 
(a,0,0,p,q,v)  distribution;  no  goodness-of-fit  test  is  attempted 
for  this  data  set  since  Arbous  and  Sichel  do  not  show  their 
grouping.  They  report  a x value  of  17.0  with  13  degrees  of 
freedom  (df)  and  P=0.20  associated  with  the  G-B-N(l,0,v)  model 
so  certainly  the  fit  is  reasonable.  A visual  comparison  of 
the  expected  cell  frequencies  for  the  two  models  indicates  a 
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close  agreement  and  so  we  would  expect  a similar  probability 
P to  obtain  for  the  fit  of  the  BNB (a  ,0 , 0 ,p  ,q ,v)  distribution. 

Although  the  fit  of  the  G-B-N(l,0,v)  model  to  the  ob- 
served data  is  reasonably  good,  the  authors  point  out  that 
12  of  18  observed  means  lie  below  the  theoretical  regression 
function  CMLE  estimates  are  used  in  (2.12)).  The  BNB  (a  ,0 , 0 ,p  ,q  ,v) 
model,  and  using  MLE  estimates,  gives  rise  to  a regression  func- 
tion for  which  only  10  of  the  18  observed  means  are  less  than 
the  predicted  values. 

We  tried  to  fit  the  BNB (a ,b ,c ,p ,q , 1)  and  the  BNB (a ,b ,0 ,p ,q ,v) 
models  in  equations  (2.36)  and  (2.37),  respectively,  to  these 
data  but  got  zero  MLE  estimates  of  b and  c in  the  first  model 
and  an  estimate  of  b equals  zero  in  the  second.  The  apparent 
lack  of  influence  of  the  parameters  b and  c will  be  noted  again 
for  the  Bates -Neyman  data.  No  attempt  was  made  to  fit  the  BNB-TR 
(0,v^,V2,Vj)  distribution  of  (2.57)  to  this  data  set  since  the 
marginal  estimates  of  the  parameter  v (column  4 of  Table  1)  are 
approximately  the  same.  Next  we  discuss  the  disease  data. 

In  their  paper  Bates  and  Neyman  present  several  data  sets 
related  to  injuries  and  diseases  suffered  during  a common  period 
of  time  by  office  and  industrial  workers.  For  the  set  associ- 
ated with  two  kinds  of  diseases  for  1286  industrial  workers  the 
fit  of  the  G-B-N(a,0,v)  distribution  is  deficient  both  in  des- 
cribing the  bivariate  data  and  the  observed  regression  function. 
That  a bnb  distribution  should  be  a candidate  model  for  the 
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data,  though,  is  suggested  by  Table  1 since  the  marginals  are 
reasonably  well  fitted  by  a univariate  negative  binomial  and 
the  sample  correlation  coefficient  is  moderate. 

In  the  following  diagram  we  show  the  empirical  and  theo- 
retical regression  functions  for  this  data;  theoretical  values 
result  by  substituting  MLE  parameter  estimates,  via  fitting 
the  joint  raw  data,  into  (2.12).  Although  a "by-eye"  fit  of 
a regression  function  to  the  data  fails  somewhat  due  to  the 
fact  that  the  observed  means  are  based  on  varying  sample  sizes, 
shown  at  the  bottom  of  the  diagram,  it  is  apparent  that  there 
is  some  nonlinearity  in  the  data. 


THEORETICAL  AND  OBSERVED  REGRESSION 
FUNCTIONS  FOR  BATES -NEYMAN  DATA 
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Table  4 shows  the  observed  and  expected  cell  frequencies 
given  by  Bates  and  Neyman  for  the  disease  data;  the  expected 
frequencies  are  based  upon  the  G-B-N  (a , 0 ,v)  distribution.  From 
the  previous  diagram  and  the  table  it  is  clear  that  this  par- 
ticular model  fails  to  adequately  describe  the  joint  data  and 
the  regression  function  so  it  is  natural  to  seek  a better  rep- 
resentation for  the  data.  To  this  end  we  fit  to  the  data  the 
new  bnb  distributions;  our  results  are  not  completely  success- 
ful in  that  the  joint  fits  are  also  inadequate  even  though  the 
X values  are  much  reduced  from  that  of  the  G -B-N (a , 0 ,v)  fit. 

That  the  empirical  regression  function  is  described  better  will 
be  illustrated. 

We  now  give  in  displays  like  Table  4 the  results  of  apply- 
ing the  new  bnb  distribution  to  these  data.  In  order  to  pro- 
vide a close  comparison  with  the  Bates-Neyman  results  we  change 
their  groupings  only  when  necessary  to  maintain  the  minimum  ex- 
pectation of  three.  Tables  5 through  8 show  the  results  and 
Figure  1 gives  some  of  the  associated  regression  functions. 

Table  5 illustrates  our  first  attempt  at  fitting  these  data; 

the  bivariate  geometric  was  chosen  since  it  leads  to  a nonlinear 

2 

regression  function  (see  2.44).  Choosing  a criterion  of  x > the 
fit  is  much  better  than  the  one  in  Table  4 but  still  inadequate. 
Figure  la  gives  the  regression  function  from  (2.44)  with  MLE 
estimates;  lb  shows  the  least  squares  fit  of  the  curve  dg+c^k*, 
which  is  the  form  of  (2.44),  to  the  data.  We  emphasize  the 
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fact  that  all  parameter  estimates  are  based  on  the  raw  data; 
that  is,  these  curves  were  not  fitted  to  the  illustrated  means. 

Tables  6 and  7 show  how  the  BNB(a,b ,0 ,p ,q,v)  in  (2.37) 
and  the  BNB (a , 0 , 0 , p ,q , v)  in  (2.26),  respectively,  describe  the 
bivariate  data.  Figures  lc  and  d give  the  corresponding  re- 
gression functions.  We  note  that  the  regression  function  for 

the  BNB(a ,b  ,0  ,p ,q ,v)  distribution  is  nonlinear.  It  is  observed 
2 

that  the  x values  are  approximately  the  same  in  Tables  5,  6 
end  7 so  none  of  the  special  cases  considered  seem  superior  in 
describing  the  joint  observations.  The  regression  functions 
are  different,  though,  and  it  appears  as  if  the  BNB (a ,b , 0 ,p ,q , v) 
regression  model  is  best,  at  least  among  the  ones  based  on  MLE 
estimates . 

A possible  reason  for  the  difficulty  in  adequately  fitting 
these  data  is  that  each  of  the  distributions  discussed  in  Tables 
4 through  7 has  a common  v parameter  associated  with  the  margin- 
als, but  the  individual  sample  values  of  v from  Table  1 are 
quite  different.  For  one  marginal  the  estimate  of  v is  0..53 
and  for  the  other,  1.69.  Thus  we  are  lead  to  apply  the  BNB-TR 
(0,v1,V2,Vj)  distribution  in  (2.57);  Table  8 displays  the  fit. 

We  see  an  improvement  in  the  x 2 value  but  still  not  enough  to 
produce  a reasonable  fit.  Figure  le  shows  the  regression  func- 
tion for  this  model.  For  reference  we  give  the  least  squares 
linear  regression  in  Figure  If. 
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Associated  with  each  bivariate  display  in  Tables  4 to  8 
is  an  implied  marginal  fit.  None  of  the  distributions  ade- 
quately describe  both  observed  marginal  distributions  although 
the  G-B -N (a , 9 ,v)  model  does  describe  adequately  (P=0.30)  the 
random  variable  labeled  respiratory  disease  and  the  BNB-TR 
(9 *^2 3)  model  gives  P=0.05  for  the  digestive  disease  mar- 
ginal and  P“0.02  for  the  other  one. 

For  the  Arbous -Sichel  data  and  for  special  cases  of  the 
general  BNB(a,b,c,p,q,v)  model  we  saw  that  the  parameters  b 
and  c were  not  needed  to  describe  those  data.  Here  for  the 
Bates -Neyman  data  we  see  from  Tables  6 and  7 that  the  impact 
of  the  parameter  b appears  minimal  in  that  the  expected  cell 
frequencies  are  about  the  same.  Ignoring  the  observation  that 
the  BNB(a,b ,0 ,p ,q , ) model  leads  to  a nonlinear  regression 
function,  whereas  the  BNB (a ,0  ,0  ,p  ,q ,\>)  distribution  gives  a 
linear  one  as  shown  in  (2.50),  we  suspect  for  parent  populations 
with  positive  correlation  that  the  latter  model  is  fairly  ro- 
bust against  alternatives  involving  nonzero  parameters  b and 
c.  Although  not  attempted  here,  this  conjecture  could  be  ex- 
amined via  simulation;  plots  of  the  probability  surface  as  a 
function  of  some  of  the  parameters  could  be  helpful  too. 

2 . 5 Summary 

In  this  part  we  discussed  the  bnb  distribution  introduced 
by  Guldberg  (1934)  and  generalized  by  Bates  and  Neyman  (1952). 
This  bnb  distribution  admits  of  positive  correlations  and  linear 
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regression  functions.  Also  we  introduced  and  derived  several 
properties  for  two  new  bnb  distributions,  one  obtained  by  con- 
volving a bivariate  geometric  distribution  given  by  Paulson 
and  Uppuluri  (1972)  , and  another  obtained  by  reducing  a certain 
trivariate  negative  binomial  distribution.  For  the  convolution 
process  a dual  bivariate  gamma  distribution  exists  (Paulson 
(1973))  and  for  it  the  duality  implies  the  exactly  analogous 
properties. 

For  the  bnb  distribution  resulting  from  a convolution, 
labeled  BNB(a,b,c,p,q,v)  , we  established  the  following  results: 

(1)  for  v integer,  the  probability  function  in  (2.32), 

(2)  for  v>0,  the  moments  to  order  two  in  (2.43  a and  b) , 

(3)  for  c=0  and  v>0,  the  probability  function  in  (2.37), 
the  probability  density  function  for  the  dual  bivariate  gamma 
distribution  in  (2.41)  and  the  nonlinearity  of  the  regression 
function  in  Figure  lc, 

(4)  for  v*l,  the  probability  function  in  (2.36)  and  the 
equation  for  the  regression  function  (nonlinear)  in  (2.44), 

(5)  for  b=c=0  and  v>0,  a method  of  moments  parameter  es- 
timation technique  in  (2.43  c-f),  the  equation  for  the  regres- 
sion function  in  (2.50),  a recursive  formula  for  the  probabili- 
ties in  (2.51),  the  likelihood  equations  for  a random  sample 
(assume  v known)  in  (2.52),  the  differential-difference  equa- 
tions (v  known)  in  (2.53),  in  (2.54)  a necessary  relationship 
for  the  parameters  in  optimizing  the  likelihood  function  and 
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the  likelihood  equations  for  the  dual  bivariate  gamma  in  (2.55) 
and  the  necessity  condition  for  its  parameters  in  (2.56). 

Contrary  to  the  Guldberg-Bates-Neyman  distribution  this 
one  gives  rise  to  positive  or  negative  correlations  and  to 
linear  or  nonlinear  regression  functions. 

For  the  bnb  distribution  which  resulted  from  a reduction 

we  derived  the  probability  function  in  (2.57),  the  moments  to 

order  two  in  (2.58)  and  the  regression  function  in  (2.59).  This 

distribution  is  basically  different  from  the  two  previously 

discussed  ones  in  that  its  marginals  have  characteristic  func- 

- v • 

tions  of  the  form  [l+0(l-elt)]  x,  i=l ,2, whereas  the  character- 
istic functions  associated  with  the  latter  distributions  are  of 
the  form  [ 1^  6 ^ (l-eit:)  ] "v  , i=l,2.  Equation  (2.1)  shows  how  the 
resulting  marginal  distributions  would  differ.  From  the  appli- 
cations of  these  distributions  to  data  we  suspect  that  bnb 
models  which  give  rise  to  marginals  with  characteristic  func- 

-V  • 

tions  of  the  form  [l+0i (l-e1*) ] 1,  i = l,2,  would  be  useful. 

We  applied  these  distributions  to  the  bivariate  data  given 
by  Arbous  and  Sichel  (1954)  on  absenteeism  among  248  shift  wor- 
kers in  two  yearly  periods  and  to  disease  data  of  two  types 
among  1286  industrial  workers. 

Arbous  and  Sichel  fitted  the  symmetric  (a=l)  Guldberg- 
Bates-Neyman  model  in  (2.6)  to  absenteeism  data  and  got  a rea- 
sonable fit  (P^O^O)  but  the  regression  function,  with  para- 
meters via  MLE,  overestimated  the  observations  in  that  12  of 
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18  means  were  below  the  computed  regressions.  We  applied 

the  BNB (a , 0 , 0 ,p ,q  , v)  distribution  to  these  joint  data  and 

the  regression  function  was  such  that  only  10  of  18  means 

were  below  it.  Arbous  and  Sichel  did  not  show  the  cell 

2 

groupings  they  adopted  for  use  in  the  x test  for  the  bivari- 
ate fit  so  we  could  not  compare  the  BNB  (a , 0 , 0 ,p  ,q  , v)  model 
to  the  Guldberg- Bates -Neyman  model.  By- eye,  the  fits  seemed 
comparable . 

Bates  and  Neyman  applied  the  Guldberg- Bates -Neyman  dis- 
tribution to  disease  data  but  the  fit  did  not  adequately  des- 
cribe the  observed  bivariate  data  or  the  observed  means  which 
clearly  suggested  a nonlinear  form  for  the  regression  func- 
tion. We  applied  several  special  cases  of  the  BNB  (a  ,b  ,c  ,p ,q , v) 
model  and  the  so  called  BNB-TR(0 ,v^ , ,Vj)  model  in  (2.57)  to 

these  data  and  never  got  a reasonable  fit  although  the  values 

2 2 
of  x were  much  reduced;  from  the  latter  model  the  x value 

was  about  one-half  of  the  value  reported  by  Bates  and  Neyman. 
The  nonlinear  regression  functions  resulting  from  the  special 
cases  of  the  BNB (a ,b ,c ,p ,q ,v)  distribution  described  much  bet- 
ter the  observed  means. 

From  our  experience  in  fitting  special  cases  of  the  BNB 
(a  ,b  ,c  ,p  ,q,v)  distribution  to  these  data  and  to  the  data  to  be 
discussed  in  the  next  part  we  suspect  that  the  parameters  b 
and  c are  relatively  unimportant  in  fitting  bivariate  data 
from  populations  with  positive  correlation.  Except  in  those 
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cases  where  the  empirical  regression  function  is  obviously 

♦ 

nonlinear  the  BNB (a , 0 , 0 ,p  ,q  , v)  case  is  probably  a fairly  ro- 
bust model.  The  recursive  formula  and  parameter  reduction 
technique  developed  here  make  it  relatively  easy  to  work  with, 


too . 


PART  3 


SOME  BIVARIATE  APPROACHES  FOR  ANALYZING  AIRCRAFT 
OPERATIONS  AND  MAINTENANCE  DATA 

3 . 1 Introduction  and  Historical  Review 

Our  objective  in  this  part  is  to  apply  bivariate  distri- 
butions to  some  problems  related  to  inventory  control  and 
maintenance  in  military  aircraft  logistics.  In  a way  to  be 
shown  we  form  bivariate  r.v.'s  related  to  these  problems  and 
illustrate  their  utility  with  actual  data.  Although  the  dis- 
cussion is  restricted  to  applying  two  of  the  aforementioned 
bnb  distributions  to  certain  data  sets,  the  techniques  are 
applicable  to  other  settings. 

The  analyses  presented  here,  for  the  most  part,  are  in 
the  context  of  fitting  distributions  to  bivariate  data  and 
then  using  these  sample  distributions  to  address  certain  prob- 
lems. Several  authors  have  postulated  areas  in  reliability 
where  certain  continuous  bivariate  distributions  can  be  ex- 
pected to  result.  See,  for  example,  Downton  (1970),  Harris 
(1968),  Hawkes  (1972),  and  Marshall  and  Olkin  (1967)  who  study 
bivariate  exponential  distributions.  Closer  to  the  techniques 
envisioned  here  are  the  works  of  Fawcett  and  Gilbert  (1966)  for 
characterizing  demand  patterns  (univariate)  for  aircraft  spare 
parts  and  Youngs,  Geisler  and  Brown  (1955)  for  predicting  de- 
mand for  aircraft  spare  parts  using  the  method  of  conditional 
probabilities.  The  account  by  Dade  (1973)  for  some  alternative 
approaches  to  maintenance  analysis  is  of  interest  also. 
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In  the  next  section  we  present  some  bivariate  approaches 
to  problems  dealing  with  discrete  data.  Section  3 describes 
how  similar  approaches  could  be  used  with  continuous  data. 

3 . 2 Some  Bivariate  Analyses  of  Discrete  Data 

The  purpose  of  this  section  is  to  illustrate  with  several 
examples  how  discrete  bivariate  r.v. 's  can  be  associated  with 
certain  aircraft  operations  and  maintenance  problems  and  sub- 
sequently give  rise  to  meaningful  analysis  techniques.  In  par- 
ticular, we  present  demand  data  for  aircraft  spare  parts  and 
use  bivariate  distributions  to  suggest  how  a particular  kind 
of  inventory  model  can  be  constructed;  additionally,  aircraft 
abort  data  are  given  and  it  is  shown  how  the  regression  func- 
tion for  a certain  bivariate  r.v.  related  to  these  data  can  be 
used  to  suggest  the  effect  of  overhaul  on  an  aircraft's  per- 
formance. Besides  the  above  applications  an  important  observa- 
tion in  its  own  right  is  that  these  data,  properly  defined, 
can  be  described  adequately  by  univariate  and  bivariate  nega- 
tive binomial  distributions. 

First  we  show  that  bnb  distributions  adequately  describe 
demand  data  for  aircraft  spare  parts  in  two  adjacent  time  per- 
iods. Table  9a  gives  for  a random  sample  of  72  aircraft  parts 
actual  demand  data  for  a four  month  period  where  we  have  formed 
bivariate  data  by  splitting  the  period  into  two  smaller  inter- 


vals; we  take  the  first  interval  to  be  the  first  three  months 
and  the  second  interval  to  be  the  fourth  month.  As  an  illustra- 
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tion  of  the  data  two  parts  were  demanded  five  times  in  the 
first  three  months  and  then  a single  time  in  the  fourth  month. 
From  Table  1 we  see  that  the  sample  correlation  coefficient  for 
the  bivariate  data  is  0.54  and  that  the  marginals  are  fitted, 
reasonably  well,  with  the  univariate  negative  binomial  distri- 
bution (illustrated  in  9b).  We  expect  then  that  one  of  the 
bnb  distributions  should  describe  these  data  and  in  Table  9a 
we  show  the  fit  of  the  Guldberg-Bates -Neyman  model  of  (2.6). 

(In  this  entire  part  procedures  for  estimating  parameters,  ap- 
2 

plying  the  x test  and  illustrating  results  are  all  the  same 
as  were  described  in  Section  2.4.)  For  this  fit  we  have  P=0.18 
indicating  that  approximately  18  of  100  fits  would  be  worse, 
assuming,  of  course,  that  the  G-B-N(ot , 6 , v)  model  is  the  under- 
lying parent  population.  That  another  bnb  distribution  describes 
these  data  is  illustrated  next. 

In  Table  10  we  show  how  the  BNB (a ,0 ,0 ,p ,q , v)  distribution 
in  equation  (2.26)  describes  the  data;  this  is  the  same  data 
set  that  was  used  before  to  illustrate  how  the  P value  associated 
with  the  test  is  affected  by  cell  groupings.  We  reiterate 
that  an  alternate  goodness-of- fit  test,  one,  perhaps,  being  in- 
dependent of  cell  groupings,  would  be  desirable  but  here  we  use 

2 

the  x test  and  agree  to  report  results  based  upon  a single  at- 
tempt at  grouping  the  cells  for  the  minimum  expectation  of  three. 

By  inspection  the  fits  illustrated  in  Tables  9 and  10  seem 
to  be  about  the  same.  The  P values  associated  with  the  marginal 
fits  of  the  BNB (a ,0 ,0 ,p ,q ,v)  distribution  to  the  univariate  ob- 
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servations  are  0.24  and  0.27  and  the  G-B-N(a,6,v)  model  gives 
rise  to  similar  values. 

That  bnb  distributions  can  be  applied  to  these  kinds  of 
data  should  be  useful  in  inventory  control  and  particularly 
to  problems  where  short  range  predictions  are  required,  such 
as  in  constructing  fly-away  kits.  In  military  air  operations 
certain  parts  are  set  aside  from  normal  operating  channels  and 
in  an  emergency  are  supposed  to  provide  enough  spares  to  last 
for  a fixed  amount  of  time,  usually  one  month.  These  parts 
make  up  a so  called  fly-away  kit.  Youngs,  et.  al.  (1955)  sug- 
gest this  application  in  their  report  but  show  no  distributions 
fitted  to  bivariate  data  as  we  do.  We  leave  this  area  and  next 
discuss  some  aircraft  abort  data. 

That  bivariate  distributions  which  admit  of  negative  corre- 
lations can  be  useful  in  applications  is  illustrated  with  the 
following  data  set.  We  show  in  Table  11  flight  aborts  (missions 
interrupted  during  flight)  for  a random  sample  of  109  aircraft 
for  two  consecutive  six  month  time  periods.  The  flight  aborts 
are  limited  to  those  caused  by  materiel  failures.  From  Table  1 
the  sample  correlation  coefficient  is  -0.16  and  we  see  that  one 
of  the  marginals  is  very  well  fitted  with  the  univariate  nega- 
tive binomial  but  the  other  fit  is  inadequate.  The  expected 
cell  frequencies  in  Table  11  are  from  the  BNB(a ,b ,c ,p ,q , 1)  dis- 
tribution and  a P value  of  0.12  results. 
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Next  we  give  other  aircraft  abort  data  and  show  an  anal- 
ysis which  suggests  the  effect  of  overhaul  on  an  aircraft's 
performance.  Aircraft  undergo  large  scale  overhaul  programs 
periodically  and  it  is  important  to  know  if  the  programs  are 
beneficial.  Traditionally,  these  programs  have  been  justified 
by  the  common  assumption  that  the  aircraft  are  restored  to  a 
better  condition;  the  following  diagram  depicts  one  way  this 
benefit  is  perceived. 


maintenance 
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Months  before  Months  after 

overhaul  overhaul 

Before  overhaul,  for  an  individual  aircraft,  we  show  an  in- 
creasing trend  to  denote  degradation  of  a variable  such  as 
number  of  failures,  number  of  aborts  or  perhaps,  unscheduled 
maintenance  manhours  and  an  improvement  immediately  after  the 
event  followed  by  degradation  again. 

Previous  attempts  at  investigating  the  effect  of  overhaul 
have  centered  on  collecting  sample  data  for  some  variable  on 
either  side  of  the  event  and  then  performing,  say,  an  analysis 
of  variance  or  a regression  analysis  to  determine  if  such 


degradation  trends  do,  in  fact,  seem  reasonable.  Certain 
studies  have  indicated  that  aircraft  are  not  in  a degraded 
condition  prior  to  overhaul  and  are  not  improved  by  the  event 
(Dade  (1973)).  Some  qualifying  statements  are  required  here-- 
certainly  an  analysis  of  this  nature  is  very  complex  in  that 
many  factors  such  as  environment,  previous  missions,  flying 
hour  history  and  total  age  of  the  aircraft  are  involved  so  evi- 
dently no  one  analysis  can  be  expected  to  be  complete  and  ex- 
haustive. These  types  of  analyses  can  be  suggestive  of  the 
effect  of  overhaul,  though,  and  in  this  respect  we  seek  to 
contribute  another  technique  and  illustrate,  with  some  new 
data,  its  use.  Although  not  attempted  here  common  experimental 
design  techniques  could  be  employed  to  control  some  of  these 
other  factors. 

We  are  primarily  interested  in  developing  techniques  ap- 
plicable to  non-normal  data  and  particularly  to  discrete  data, 
such  as  aborts  which  are  typically  small.  From  here  on  we 
analyze  total  aborts,  which  are  mission  interruptions  discovered 
during  pre-flight  or  in-flight  operations.  As  before  we  study 
only  those  aborts  caused  by  materiel  failures. 

To  examine  the  effect  of  overhaul  we  define  two  bivariate 
r.v.'s  related  to  the  periods  shown  in  the  following  diagram. 
Periods  1 and  2 (common  to  all  aircraft)  are  two  adjacent  six 
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month  periods  where  no  overhaul  was  performed  and  periods  3 
and  4 are  the  six  month  periods  immediately  before  and  after 
an  overhaul,  respectively.  The  period  lengths  were  chosen 
arbitrarily.  If  is  the  number  of  total  aborts  per  aircraft 
in  period  i,  i=l,2,3,4,  then  we  wish  to  compare  the  probability 
distribution  of  the  pair  (Aj,A^)  where  an  intervening  overhaul 
was  performed  to  the  probability  distribution  of  the  pair  (A^, 
where  no  intervening  overhaul  occurred.  Under  a null  hypo- 
thesis of  the  overhaul  being  ineffective  in  changing  an  air- 
craft's performance,  as  measured  by  aborts,  these  two  distri- 
butions should  be  the  same.  The  regression  function  is  a 
descriptor  of  a bivariate  distribution  (Kendall  and  Stuart 
(1973))  and  so  one  simple  way  to  compare  these  two  distributions 
would  be  to  examine  the  two  regression  functions  associated  with 
sample  data.  All  other  factors  being  equal  any  two  aircraft 
with  the  same  number  of  aborts  in  periods  1 and  3 should  have 
the  same  number  of  aborts,  on  the  average,  for  periods  2 and  4 
if  overhaul  is  ineffective  and  so  any  difference  in  the  corre- 
sponding regression  functions  should  be  suggestive  of  the  effect 
of  overhaul. 

This  type  of  analysis  presupposes  the  existence  of  a suit- 
able bivariate  distribution  which  will  adequately  describe  ob- 
served data.  Next  we  present  sample  data  related  to  the  bivari- 
ate r.v.'s  (Aj , A2)  and  (Aj,A4)  and  fit  them  with  the  BNB 
(a,0 ,0 ,p,q,v)  distribution  in  equation  (2.26).  This  distribu- 
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tion  was  chosen  arbitrarily  and  we  suspect  that  other  bnb 
distributions  could  be  used  equally  well.  Sample  values  of 
Aj  will  be  denoted  by  aj  , j=l,2,3,4. 

Table  12  gives  the  number  of  total  aborts  in  periods  1 
and  2 for  203  aircraft.  To  obtain  this  sample  we  considered 
an  entire  population  of  aircraft  of  a particular  type  (about 
500)  and  excluded  those  aircraft  with  an  overhaul  during  the 
time  periods  specified  as  1 and  2 and  also  those  aircraft  with- 
out 12  full  months  of  reported  abort  data  during  that  time  -- 
203  aircraft  resulted.  As  an  illustration  of  the  data  three 
aircraft  had  ten  total  aborts  in  period  1 followed  by  six  total 
aborts  in  period  2. 

For  these  data  we  can  use  the  ordinary  sign  test  (Gibbons 
(1971))  for  a bivariate  r.v.  to  conclude  that  where  no  inter- 
vening overhaul  is  involved  aircraft  incur  the  same  number  of 
total  aborts  in  two  adjacent  six  month  periods  (here  we  are 
actually  testing  the  hypothesis  that  the  median  of  the  r.v. 
(A1-A2)  is  zero  versus  a two-sided  alternative;  a normal  ap- 
proximation to  the  binomial  used  in  this  nonparametric  test 
gives  a sample  value  of  |z|=1.25). 

Table  1 shows  certain  descriptive  statistics  for  these 
data  and  it  is  apparent  that  the  univariate  negative  binomial 
distribution  can  be  used  to  describe  the  marginal  observations. 
We  wish  to  describe  further  the  bivariate  data,  though,  and 
so  we  attempt  to  fit  the  joint  observations  with  a bivariate 
negative  binomial  distribution. 


so 


Table  13  shows  the  expected  cell  frequencies  which  result 

by  applying  the  BNB (a ,0 , 0 ,p ,q , v)  distribution  of  (2.26)  to  these 

data  and  Table  14  gives  the  observed  and  expected  frequencies 

2 

together  and  the  x values.  For  convenience,  we  sum  the  obser- 
vations in  any  one  "super-cell"  as  is  our  custom  for  the  expec- 

2 

tations.  We  conclude  that  the  bivariate  fit,  as  measured  by  x » 
is  good  (P=0.65).  The  P values  for  the  associated  marginal  fits 
are  0.04  for  and  0.56  for  From  equation  (2.50)  with  MLE 

parameter  estimates  we  see  that  the  regression  function,  or  the 
mean  number  of  aborts  for  period  2 given  the  observed  number  of 
aborts  for  period  1,  is  E [A2 [ a^] =4 . 53  + 0.26  a^.  It  can  be 
shown  that  the  least  squares  regression  function  for  these  data 
is  E [A2 | aj] =4. 41  + 0.28  a^.  Next  we  show  the  same  analysis  for 
sample  data  from  periods  3 and  4;  that  is,  before  and  after  an 
ove  rhaul . 

Tables  15,  16  and  17  give  observations,  expectations  and 
2 

the  x test  for  periods  3 and  4 for  387  aircraft.  These  are 
the  same  type  of  aircraft  as  before  and  we  have  an  aircraft 
being  included  in  this  sample  if  it  has  six  full  months  of  re- 
ported abort  data  on  adjacent  sides  of  a common  overhaul  event. 
We  point  out  that  periods  3 and  4 may  be  separated  by  two  or 
three  months  which  is  usually  the  length  of  an  overhaul  for 
these  aircraft. 

If  we  apply  the  sign  test  to  the  bivariate  data  for  (Aj, 

A4)  shown  in  Table  15,  as  was  done  for  the  observed  data  for 
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(A^.A-,),  there  results  a sample  value  of  |z|=  4.51  and  we  can 
conclude  that  values  of  A^  greater  than  A^  are  more  likely  in 
these  data;  that  is,  total  aborts  per  aircraft  are  generally 
greater  after  overhaul. 

Column  5 of  Table  1 shows  how  the  univariate  negative  bi- 
nomial distribution  fits  the  marginals.  For  the  joint  data, 
the  P value  (0.12)  associated  with  the  fit  of  the  BNB(a,0,0, 
p,q,v)  distribution  is  not  as  high  as  for  the  previous  data 
but  we  assume  that  it  is  acceptable.  The  implied  marginal  fits 
give  rise  to  P values  of  0.43  for  A^  and  0.05  for  A^.  As  a 
descriptor  of  these  bivariate  data,  the  regression  function 
(using  MLE  estimates)  is  E [A4 | a^] =7 . 13  + 0.26  a^.  The  least 
squares  regression  function  is  E [A^ | a^] = 7 . 32  + 0.23  a^.  Figure 
2 shows  on  one  graph  the  estimated  regression  functions  (via 
MLE)  for  these  two  data  sets.  We  view  the  regression  functions 
being  useful  in  the  following  way.  If  two  aircraft  have  the 
same  number  of  aborts  in  periods  1 and  3,  say,  ten,  then  the 
aircraft  which  has  an  overhaul  will  have,  on  the  average,  about 
nine  and  one-half  aborts  in  the  next  six  months  and  the  air- 
craft that  does  not  get  an  overhaul  will  have  about  seven,  again 
on  the  average.  An  important  advantage  in  being  able  to  fit 
the  data  as  we  have  illustrated  here  is  that  a confidence  in- 
terval could  be  placed  on  the  predicted  number  of  aborts  by 
using  the  estimated  bivariate  probability  function;  thus,  we 
avoid  the  normality  assumption  that  is  traditionally  involved 
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with  placing  confidence  limits  on  the  predictions  (Draper  and 
Smi th  (1966) ) . 

Apart  from  the  development  described  here  using  bnb  dis- 
tributions and  if  we  were  willing  to  invoke  the  necessary  nor- 
mality assumptions,  an  alternative  approach  would  be  to  use 
the  least  square  regression  functions  in  the  above  way.  The 
results  for  these  data  would  be  about  the  same. 

Of  interest  is  another  bivariate  r.v.  which  is  formed  by 
taking  and  A^  where  A^  is  a r.v.  representing  the  total 
aborts  in  the  six  months  immediately  after  the  fourth  period. 
Thus,  we  have  a new  bivariate  r.v.  associated  with  the  aircraft's 
performance  right  after  overhaul  and  we  are  interested  in  how 
it  compares  with  (A3,A4),  the  r.v.  representing  total  aborts 
on  either  side  of  overhaul.  Specifically,  we  are  interested 
in  whether  or  not  aborts  decrease  again;  if  the  overhaul  is 
responsible  for  the  observed  increase  in  period  4 perhaps  the 
degradation  is  similar  to  that  usually  experienced  in  a new 
item’s  performance  and  which  gradually  declines  to  a lower 
level . 

Table  18  shows  sample  data  for  (A4,A5)  from  428  aircraft 
of  the  same  type  previously  considered;  forming  the  differ- 
ences (a^-ajj)^,  i*l,2, ...  ,428,  and  using  the  sign  test  for  an 
hypothesis  of  zero  median  difference  for  the  r.v.  A4*A,-,  versus 
a two-sided  alternative,  leads  to  a sample  value  of  z=0.74  and 
so  we  take  the  hypothesis  to  be  true  for  these  data.  Evidently 
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then,  the  number  of  total  aborts  is  about  the  same  for  periods 
4 and  5 for  an  individual  aircraft. 

We  point  out  the  rather  large  values  of  P,  0.75  for  A4 
and  0.56  for  A^  , which  result  from  fitting  the  univariate  nega- 
tive binomial  distribution  to  the  marginal  observations  (see 
column  5 of  Table  1) . Although  not  illustrated,  we  applied  to 
these  sample  data  the  same  bnb  distribution  that  was  used  in 
the  two  previous  instances;  the  MLE  parameter  estimates  which 
resulted  are  a=0.7485,  p=0.7728,  q=0.7642  and  v=2.0285.  Assoc- 
iated with  the  bivariate  fit  is  P=0.27  and  for  the  marginals, 
P=0.66  for  A4  and  P=0.76  for  A,-.  We  accept  the  fit.  Based  on 
MLE  estimates,  the  regression  function  is  E [A5 | a4] =6 . 97  + 0.20  a4 
and  from  least  squares,  E [A^ | a4] =6 . 84  + 0.21  a4 ; the  dashed  line 
in  Figure  2 shows  the  former.  The  upper  two  regression  func- 
tions in  the  figure  are  probably  within  sampling  error  of  one 
another;  no  attempt  is  made  to  test  for  true  differences. 

Although  not  done  here  it  would  be  of  interest  to  make  a 
similar  comparison  before  overhaul;  that  is,  take  a six  month 
period  preceding  period  3 and  form  a bivariate  r.v.  for  the 
number  of  total  aborts  in  the  new  period  and  period  3 and  then 
compare  that  r.v.  to  (A3,A4).  Another  analysis  planned  for 
the  future  is  to  compare  ^^^2)  to  (Aj,A4)  where  just  those 
aircraft  with  overhauls  during  periods  1 and  2 are  selected 
for  consideration  in  periods  3 and  4.  In  this  way  we  have  two 
samples  during  the  same  calendar  time  (approximately)  and  any 


possible  contributing  causes  from  that  factor  would  be  some- 
what controlled;  Dade  points  out  that  changes  in  management 
policy  or  hardware  and  tactics  connected  to  calendar  time  have 
often  made  maintenance  data  analyses  difficult. 

For  a large  group  of  aircraft  of  a particular  type  our 
results  show  how  these  aircraft  performed  after  overhaul.  The 
results  are  only  suggestive;  with  only  the  surface  analysis  des- 
cribed here,  we  are  not  prepared  to  argue  that  the  observations 
represent  true  degradation  or  that  they  were  caused  by  the  over- 
haul. We  do  think  the  approach  is  worthy  of  consideration 
though  and  coupled  with  a comprehensive  experimental  design 
could  provide  unbiased  results. 

Other  definitions  of  the  periods  might  lead  to  more  re- 
vealing results;  some  alternatives  are  fixed  flying  hour  inter- 
vals, a fixed  number  of  sorties,  or  shorter  monthly  intervals. 
Here  the  total  aborts  were  taken  for  the  whole  aircraft;  per- 
haps total  aborts  for  a particular  subsystem  would  be  better. 

In  the  next  section  we  discuss  briefly  similar  applications 
for  problems  dealing  with  continuous  data. 

3.3  Continuous  Bivariate  Distributions  Applied  to  Aircraft 

Failure  Data 

In  this  section  we  describe  how  bivariate  probability  den- 
sity functions  associated  with  certain  aircraft  failure  data 
could  be  useful  in  managing  aircraft  operations  and  maintenance 
programs.  Although  the  discussion  is  centered  on  using  bivari- 
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ate  gamma  distributions  the  techniques  are  applicable  to 
other  bivariate  densities.  First  we  describe  the  basic  ran- 
dom variable  of  interest  and  then  show  how  meaningful  bivari- 
ate models  can  be  constructed. 

For  an  aircraft  we  take  as  a r.v.  X the  operating  time 
between  failures  on  any  component  for  which  the  univariate 
gamma  distribution  is  expected  to  adequately  describe  the  ob- 
served failure  times.  An  examination  of  historical  data  would 
indicate  candidate  components.  (We  use  the  term  component  to 
represent  an  individual  part  or  a subsystem  of  parts.)  For 
applications  to  be  defined  later  we  also  require  on  each  air- 
craft the  delivery  date  and  dates  of  overhauls.  These  over- 
hauls are  on  the  entire  aircraft  and  not  just  on  the  component 
itself. 

On  the  following  time  axis  we  illustrate,  for  a particu- 
lar aircraft  and  component,  the  random  variable  and  events 
described  above. 

Ove rhaul 

Delivery  x 2 i-1  i i+1  -2  -1  +1  +2 

r • 1 * +-■...  1 1 1—  — t 1 1 1 time 

X1  x2  Xi  Xi  + 1 X B XA 

Failures  occur  at  epochs  1,2,...,  and  the  observed  operating 

times  are  Xj,X2 Here  the  observation  x^  is  taken  to  be 

a realization  of  which  is  the  random  variable  representing 
the  operating  time  between  the  (i-l)st  and  ith  failures.  An 
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aircraft  overhaul  occurs  at  the  indicated  point  and  the 
previous  failures  are  labeled  -1  and  -2;  the  intervening 
erating  time  is  denoted  by  Xg . In  a like  manner  the  nex 
failures  after  overhaul  are  labeled  +1  and  +2  and  the  in 
vening  operating  time  is  x^.  Next  we  describe  some  biva 
models  and  possible  uses. 

For  some  components  it  is  reasonable  to  expect  that 
operating  times  and  X^  + j are  related  (particularly  fo 
ponents  that  are  repaired  and  replaced  on  the  same  aircr 
and  so  we  assume  the  existence  of  a bivariate  distributi 
describe  this  dependence.  Given  bivariate  data  (xi>xi+i 
n=l,2,...,N,  we  could  determine  if  the  sample  values  x^ 
x^+1  are  correlated  and  if,  in  fact,  the  observed  data  c 
be  fitted  adequately  with  any  of  the  existing  bivariate 
distributions.  To  describe  some  possible  uses  of  bivari 
probability  models  in  this  context  we  assume  the  data  co 
be  fitted  adequately. 

If  the  random  variables  and  X^  + ^ can  be  describe 
a bivariate  distribution  and  if,  for  a particular  aircra 
an  observation  is  given,  then  the  conditional  distrib 
of  X^  + jJx^  could  be  used  to  predict  the  time  of  the  next 
on  this  aircraft,  at  least  due  to  the  component  being  co 
The  regression  function,  that  is,  the  mean  of  the  condit 
distribution  could  be  used  to  predict  the  mean  time  of  t 
failure.  This  would  be  an  improvement  over  a prediction 
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equal  to  the  mean  time  between  failure  (MTBF)  which  is  typi- 
cally used  and  is  computed  as  if  the  random  variables  and 
Xi+1  are  independent  and  identically  distributed. 

Aircraft  components  generally  are  assumed  to  exhibit 
failure  rates  which  vary  with  time  according  to  the  classic 
bath  tub  curve.  Thus,  failures  occur  more  rapidly  at  first 
then  decrease  to  a somewhat  constant  level  and  finally  increase 
in  frequency.  Failures  in  these  three  periods  commonly  are 
called  initial,  chance,  and  wear-out  failures,  respectively 
(Mann,  et.  al.  (1974)).  A comparison  of  the  bivariate  den- 
sities associated  with  the  sample  operating  times  on  either 
side  of  several  failures  could  suggest  a pattern  of  component 
aging. 

One  of  the  prime  considerations  in  simulation  studies  in- 
volving a composite  of  aircraft  operations  and  support  func- 
tions is  in  generating  realistic  component  failures.  Certainly, 
if  times  to  failure  are  dependent  then  a mechanism  which  allows 
for  pairs  of  observations  with  the  proper  correlation  would  be 
an  improvement  over  a model  which  ignored  the  dependence. 

A major  objective  of  aircraft  overhaul  is  to  restore  the 
aircraft  to  a more  reliable  condition.  As  before  we  are  inter- 
ested in  determining  if  the  overhaul  does  improve  performance. 

If  we  assume  the  existence  of  a bivariate  model  to  describe  the 
dependence  between  and  X^  + 1 where  no  overhaul  is  involved, 
then  we  can  investigate  the  effect  of  overhaul  by  comparing  a 
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sample  distribution 
bution  which  results 
from  a random  sample 
These  are  but  a 
tions  which  could  be 


associated  with  (X^,X^+^)  to  the  distri- 
by  fitting  the  observed  data  (xg,x^) 
of  aircraft  with  an  overhaul, 
few  of  the  bivariate  r.v.’s  and  applica- 
described  for  continuous  data. 


3 . 4 Summary 

In  this  part  we  applied  discrete  bivariate  probability 
distributions  to  actual  aircraft  operations  and  maintenance 
data  and  showed  how  these  distributions  could  be  used  in  prob- 
lems related  to  inventory  control  and  aircraft  overhaul.  Addi- 
tionally, we  described  how  continuous  bivariate  distributions 
could  be  used  to  analyze  certain  continuous  failure  data.  All 
of  the  applications  presented  were  for  self-pairing  type  situa- 
tions; without  doubt,  bivariate  applications  exist  also  in 
situations  involving  a dependency  between  two  separate  items. 
The  latter  is  a more  traditional  approach  for  continuous  bi- 
variate r.v. 's . 

We  presented  several  new  data  sets  related  to  demand  for 

aircraft  spare  parts  and  materiel  failure  induced  aborts.  For 

the  demand  and  abort  data  we  investigated  the  marginal  and 

joint  data  and  showed  how  the  univariate  and  bivariate  negative 

binomial  distributions  fitted  the  observations.  The  P values 

? 

associated  with  the  x goodness-of- f it  test  ranged  between  0.01 
and  0.85  with  the  average  being  about  0.36.  One  of  these  data 


59 


sets  had  a negative  sample  correlation  coefficient  which  em- 
phasizes the  need  for  bnb  distributions  that  can  describe  such 
dependencies  (in  Part  2 we  discussed  such  a distribution) . 

As  a possible  way  to  investigate  the  effect  of  a noted 
event  on  an  item's  performance  we  illustrated,  with  aircraft 
abort  data  and  the  overhaul  event,  how  bivariate  distributions 
•could  be  used.  Our  methodology  involved  comparing  two  bivari- 
ate distributions,  one  defined  for  r.v.'s  on  either  side  of 
the  event  and  the  other  defined  for  a similar  r.v.  not  separated 
by  the  event  of  interest.  We  used  the  sample  regression  func- 
tions to  compare  the  distributions. 

Extensions  to  multivariate  (beyond  two)  settings  are  ob- 
vious. Certainly  the  difficulty  of  obtaining  parameter  esti- 
mates is  compounded  though. 

In  the  next  part  we  apply  a bivariate  exponential  distri- 
bution to  some  correlated  queueing  systems. 


PART  4 


CORRELATED  QUEUEING  SYSTEMS 

4 . 1 Introduction  and  Historical  Review 

Using  simulation  techniques  Paulson  and  Beswick  (1973) 
showed  the  effect  of  dependent  exponential  service  times  on 
queues  in  series.  In  this  part  we  review  their  work  and  pre- 
sent a set  of  recursive  formulae  useful  in  simulating  the 
queueing  process.  We  use  spectral  analytic  techniques  to  show 
that  the  effect  is  indeed  statistically  significant.  Also,  we 
investigate  the  effect  of  correlated  interarrival  and  service 
processes  on  single  server,  single  stage  queues. 

First,  we  describe  two  physical  settings  where  tandem 
queues  with  dependent  service  times  can  be  expected  to  arise. 

In  a paper  mill,  large  rolls  of  paper  typically  pass 
through  an  inspection  or  winding  operation  prior  to  being  cut 
into  smaller  rolls.  A poor  quality  roll  takes  a relatively 
longer  time  in  the  inspection  process  because  defective  sec- 
tions must  be  removed  and  splices  made.  When  this  same  roll 
reaches  the  final  cutting  stage  it  must  be  processed  more  slowly 
to  avoid  breaking  the  splices  and  to  repair  them  when  they  do 
break.  Hence  process  times  at  the  two  stages  tend  to  be  cor- 
related; indeed,  it  is  conceptually  possible  that  they  be  highly 
correlated.  The  process  times  at  the  two  stages  on  any  two  dif- 
ferent rolls  would  generally  be  independently  distributed.  In 
the  current  context  considerable  interest  would  be  centered  on 
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the  effect,  if  any,  produced  by  non- independence  of  process 
times  at  different  stages. 

Jackson  (1954)  in  discussing  queueing  systems  with  phase 
type  service  pointed  out  that  a typical  sequence  of  events  in 
the  overhaul  of  an  aircraft  engine  consist  of  stripping,  de- 
tailed examination,  repairs,  assembly  and  testing.  Generally, 
an  engine  with  a large  number  of  maintenance  requirements  can 
be  expected  to  spend  more  time  in  each  of  the  latter  four  phases 
and  so  the  possible  effect  of  correlated  service  times  on 
throughput  time  would  be  of  interest.  It  is  not  difficult  to 
envision  a host  of  other  situations  involving  queues  in  series 
in  which  the  service  times  at  the  various  stages  for  a given 
customer  are  correlated. 

A large  proportion  of  the  literature  concerning  tandem 
queues  has  centered  on  Poisson  arrival  processes,  exponential 
service  times,  and  steady  state  solutions.  The  assumption  of 
independence  of  service  times  is  intricately  interwoven  into 
the  fabric  of  the  traditional  birth-death  equation  approach  to 
finding  a transient  and  steady  state  solution  to  the  tandem 
queueing  phenomenon.  We  shall  remain  within  this  same  framework 
with  the  exception  that  we  shall  drop  the  heretofore  universal 
(but  tacit!)  assumption  of  mutual  independence  of  all  exponen- 
tial service  times.  An  obvious  approach  is  to  use  a multivar- 
iate exponential  distribution  with  non-zero  correlations  in 
place  of  the  usual  independent  exponential  service  times.  In 
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our  situation  it  is  not  clear  that  the  birth-death  equation 
approach  can  be  modified  to  incorporate  dependent  service 
times.  Moreover,  any  such  formulation  would  very  likely  be 
analytically  intractable.  The  problem  is  however  amenable  to 
a simulation  approach  and  it  is  in  this  way  that  we  assess  the 
effect  of  departures  from  independence  of  service  times  on 
steady  state  system  performance. 

The  bulk  of  this  part  is  concerned  with  two  stage  queues 
in  series  since  our  main  concern  is  showing  that  a substantial 
effect  on  system  performance  is  indeed  induced  by  correlated 
service  times. 

For  the  unusual  single  server,  single  stage  queue  where 
the  length  of  a customer's  service  time  is  determined,  with 
probability  one,  by  the  length  of  the  interarrival  interval 
separating  himself  and  his  predecessor,  Conolly  (1968)  gives 
the  waiting  time  distribution  and  its  moments.  It  is  also 
shown  that  this  pattern  of  server  behavior  results  in  a dras- 
tic reduction  of  the  mean  and  variance  of  the  waiting  time  as 
compared  with  the  conventional  M/M/1  system.  Conolly  refers 
to  this  type  of  system  as  self-regulating  and  other  results 
are  given  by  he  and  Hadidi  (1969,  1974).  We  study,  via  simu- 
lation, this  type  of  system  where  the  dependence  between  the 
service  time  and  interarrival  interval  is  assumed  to  be  proba- 
bilistic according  to  a bivariate  exponential  distribution. 
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Other  investigations  involving  dependency  between  the 
service  and  interarrival  processes  have  typically  allowed  for 
a dependence  between  a service  time  and  an  arrival  one  inter- 
val later  than  the  one  we  consider;  See,  for  example  Lindley 
(1952)  and  John  (1963).  Next  we  discuss  tandem  queues. 

4 . 2 The  Effect  of  Correlated  Exponential  Service  Times  on 
Single  Server  Tandem  Queues 

4.2.1  The  Tandem  Queueing  System  and  Recursive  Formulae 
Consider  the  tandem  queueing  process  depicted  in 
the  following  diagram.  Customers  from  an  infinite  popu- 
lation arrive  at  a two  stage  system  according  to  a Poisson 
process  with  mean  rate  X which  we  shall,  without  loss  of 
generality,  take  to  be  unity.  An  unlimited  queue  is  al- 
ways allowed  before  the  first  stage  but  before  the  second 
stage  the  queue  length  may  be  either  restricted  or  unlimited. 
A single  server  is  allowed  at  each  stage;  the  service  dis- 
cipline is  first-come,  first-served. 


CUSTOMERS  ARRIVE  IN 
ACCORDANCE  WITH  A POISSON 
PROCESS  WITH  INTENSITY  X 

cn+  2 * cn+ 1 ,cn ' cn- 1 
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CUSTOMER  SERVICE  TIMES  ARE 
GOVERNED  BY  A BIVARIATE 
EXPONENTIAL  DISTRIBUTION  WITH 
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The  system  performance  measure  is  taken  to  be  mean 
waiting  time  per  customer  and  in  this  section  we  develop 
a set  of  formulae  to  recursively  compute  the  waiting  time 
per  customer.  We  use  the  recursive  formulae  for  the  un- 
limited interstage  storage  case  in  order  to  demonstrate 
precisely  how  the  queueing  system  consisting  of  two  stages 
in  series  with  dependent  service  times  is  related  to  a 
single  server  system  with  interdependent  arrival  and  ser- 
vice processes  as  discussed  by  Bhat  (1969).  An  interpre- 
tation by  Conolly  (1968)  for  a special  type  of  this  latter 
interdependence  is  shown  to  be  helpful  in  suggesting  why 
mean  waiting  time  is  affected  by  correlated  service  times. 

Denote  by  (T  1 , T _)  the  times  between  arrival 
epochs  of  customers  c . and  c at  the  first  and  second 
stages  and  let  cn  experience  the  service  times  (Sn  Sn  2) 

at  each  stage,  n=l,2 The  sequences  of  interarrival 

times  (T  ,,  T 7)  and  the  (S  , , S_  7)  for  different  cus - 
tomers  are  both  assumed  to  be  mutually  independent  and  in- 
dependent of  each  other. 

We  take  (w  ^ , w*-2^)  to  be  the  waiting  times,  ex- 
n n 

eluding  service,  and  (W^1^ , W^2^)  to  be  the  total  waiting 
times,  of  customer  cn , at  the  respective  stages,  n=l,2,...  . 
We  illustrate  these  definitions  with  an  arbitrary  combina- 
tion of  arrival  and  service  times  in  the  following  diagram. 
The  illustration  is  for  two  queues  in  series  with  unlimited 
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interstage  storage;  diagrams  like  this  are  useful  in 
developing  the  recursive  formulae  for  the  different  cases 
to  be  presented. 
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Case  A.  Two  stage  queues  in  series,  unlimited  inter- 
stage storage 

Customer  cn+1's  total  waiting  time  at  the  first 
stage  and  interarrival  time  at  the  second  stage  are  given 
by 
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“nil  '< 


n+1,1’ 

if 

T .>Wfl) 

n+1,1—  n 

■(  1) -t  +S 

n n+1,1  n+1,1’ 

if 

T 

n+1,1  n 

(4.1) 
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n+1,2 
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The  condition  in  (4.1)  and  (4.2)  that 

T simplv  means  that  c ..  arrives  at  server 

n+1,1— n r n+1 

one  after  cn  has  departed,  and  likewise  Tn  + ^ ^ 

means  c arrives  before  c leaves. 
n+1  n 


W 


(2) 

n+1 
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to  (4.1), 

cn+l's 

wait 

ing  time  at 

stage  is 
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if 

T , 

7>w(2) 

n+1,2’ 

n+1 

, 2—  n 

w(2)_t 

if 

T , 

7<wt2). 

n n + 1 

,2  Sn+1,2’ 

n + 1 

,2  n 

(4.3) 


The  above  diagram  illustrates  (4.1),  (4.2), 
and  (4.3)  for  Tn+J  1 > w£l)  and  Tn+J  2 < W^2)  . Similar 
diagrams  result  for  the  remaining  conditions. 

In  an  obvious  way,  we  can  use  these  relation- 
ships to  build  up  a set  of  recursive  formulae  for  any 
number  of  stages  in  series  where  the  interstage  storage 
between  stages  is  unlimited  (See  Appendix  B) . 
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Since  each  customer  must  proceed  through  both 
stages  the  output  of  the  first  stage  becomes  the  input 
of  the  second  stage  and  therefore  we  have,  in  steady 
state,  that  the  time  interval  between  arrivals  at  the 
second  stage  satisfies  a Poisson  process  with  the  same 
interarrival  intensity  parameter  X as  the  input  distri- 
bution (Burke  ( 1956)).  Unlike  the  first  stage,  however, 
cn's  service  time  at  the  second  stage  is  correlated  with 
the  interarrival  time  there.  In  the  above  diagram,  this 
corresponds  to  a correlation  between  and  T .. 

This  result  is  apparent  from  (4.2)  since  Sn+1  ^ and 
Sn+j  , are  dependent  by  assumption. 

If  S . and  S ..  are  independent  as  is 
usually  assumed  for  two  stage  series  systems  then  each 
stage,  in  steady  state,  can  be  analyzed  independently  and 

since  T J.1  , and  S„.,  are  independent,  as  are  T ..  . 

n+1,2  n+1,2  r n+l,l 

and  S . , , the  regular  M/M/1  results  obtain  for  each 

n ' x f x 

stage . 


Bhat  (1969)  describes  five  different  classes  of 
single  server  first-come  first-served,  systems  with  Pois- 
son input  and  exponential  service  times  which  result  from 
relaxing  some  of  the  assumptions  of  independence  which 
are  typically  assumed.  These  classes  represent  more  rea- 
listic operating  systems  than  those  with  assumptions  of 
independence;  Bhat  further  points  out  that  more  work  needs 


I 

J 1 


68 

to  be  done  on  these  problems  than  the  limited  amount  re- 
ported at  that  time.  One  of  these  classes  is  for  systems 
with  interdependent  arrival  and  service  processes  as  is 
the  case  here  for  T - and  S 

Conolly  (1968)  and  Conolly  and  Hadidi  (1969,  1974) 
have  studied  a dependent  structure  somewhat  similar  to 
this  wherein  the  ratio  of  service  time  to  interarrival 
time  is  constant  for  all  n;  they  give  transient  as  well 
as  steady  state  results  for  the  system.  Conolly  showed 
numerically  that  this  pattern  of  server  behavior  results 
in  a drastic  reduction  in  the  mean  and  variance  of  the 
waiting  time  as  compared  with  a conventional  M/M/1  queue. 

It  was  noted  by  Conolly  that  this  kind  of  server  behavior 
is  to  be  expected  from  a well  regulated  service  facility 
where  the  server  adjusts  the  service  time  of  a customer 
according  to  that  customer's  interarrival  time,  which  the 
server  observes  without  error.  In  this  way,  a long  inter- 
val gives  rise  to  a long  service  time,  and  short  intervals 
corresponding  to  a succession  of  rapid  arrivals,  are  fol- 
lowed by  correspondingly  short  service  times.  This  regu- 
lated behavior  therefore  prevents  a long  queue  from  forming 
and  cuts  down  on  the  mean  and  variance  of  the  waiting  time 
in  the  system. 

Returning  to  the  two  stages  in  series  problem 
under  study  we  see  that  this  system,  via  equation  (4.2), 
can  be  viewed  as  a type  of  self- regulated  system  since 
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S A , , and  T are  related,  although  not  in  the  deter- 

ministic  way  assumed  by  Conolly.  It  will  be  demonstrated 

later  that  our  tvpe  of  stochastic  dependence  between  S , 

' 1 n+ 1 , Z 

and  T , 7 gives  rise  to  results  which  are  consistent  with 
Conolly's.  This  artificial  way  of  viewing  the  system  as  a 
se 1 f - regulating  device  is  employed  solely  to  make  the  ef- 
fects seem  more  reasonable  and  in  no  way  influence  the 
resul ts . 

Case  B.  Two  stage  queues  in  series,  no  interstage  storage. 

For  this  case,  c 's  total  waiting  time  at  the 
n ° 

fZl 

second  stage,  J , is  always  equal  to  the  Sn  ^ so  the 
only  quantity  of  interest  here  is  . Since  there  is 

restricted  (zero)  interstage  storage,  the  phenomenon  of 
blocking  occurs  and  so  the  waiting  time  computation  is  a 
bit  more  complicated  than  in  Case  A. 

Blocking  of  stage  one  occurs  when  a customer  who 
has  been  served  there  is  denied  entry  into  the  second 
stage  because  no  space  remains  in  the  queueing  area  for 
stage  two.  The  customer  therefore  stays  in  the  first 
stage  and  prevents  that  server  from  accepting  a waiting 
customer  for  service.  In  effect,  the  first  server's  util- 
ization is  diminished  (Saaty  (1961)). 

The  total  waiting  time  for  cn+1  at  the  first 
server  is  given  by  one  of  four  relationships  depending 
upon  the  algebraic  sign  of  - Tn+1  j --  that  is,  upon 


I 


l 


whether  or  not  cn+^  arrives  at  stage  one  before  o 

F°r  Vi,i  * wD- 


cn  leaves. 


,(1) 

n+1 


= < 


W 


CD 

n+1 


W 


r (1 ) _ 

n 

Tn  + 1 ; 

, 1+Sn+1 , 1 ’ 

if 

Sn+1 

,liS *: 

r(D_ 

n 

Tn+ 1 , 

i +s 
.1  n 

,2’ 

if 

Sn+ 1 

,1<S. 

T 

> 

w(D 

n+1 

,1  - 

n 

9 

n+1, 

1» 

if 

Sn+ 1 

(D_ 

n 

T 

n+1, 

,+S 
, 1 n 

2 9 
9 ** 

if 

Sn  + 1 

,1<W, 

(1)  T 

n n 

Cl)  _T 

n n 

The  following  diagram  illustrates  (4.4) 
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Similarly,  the  other  conditions  can  be  verified. 
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Case  C.  Two  stage  queues  in  series,  interstage  storage 
capacity  equals  one. 

As  in  the  previous  case  blocking  can  occur  at 

the  first  stage  but  here  a customer's  total  waiting  time 

at  stage  two  can  exceed  the  service  time  since  interstage 

storage  is  permitted  on  a restricted  basis. 

The  total  waiting  time  of  c at  server  one  de- 
° n+1 

pends  on  the  algebraic  sign  of  - Tr+^ 
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and  c ,, 's  interarrival  time  at  server  two  is 
n+1 
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and 
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Next  the  total  waiting  time  of  cn+1  at  stage 
two,  W^2j,  is  computed  by  using  (4.3)  in  Case  A with 
T . 7 as  defined  in  (4.7)  or  (4.9). 

The  following  diagram  is  descriptive  of  (4.8) 
and  (4.9)  where  T j^1’  and  fl^15  ‘V, _ -S„ _2 
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4.2.2  A Bivariate  Exponential  Distribution 

There  are  a number  of  bivariate  exponential  dis- 
tributions which  could  be  used  to  describe  the  dependence 
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assumed  between  and  S?  (we  drop  the  subscripts  n)  . 

We  choose  to  use  a special  case  of  the  bivariate  gamma 
distribution  discussed  by  Wicksell  (1933)  and  Kibble 
(1941),  and  more  generally  by  Krishnamoorthy  and  Partha- 
sarty  (1951)  and  Paulson  (1973).  The  functional  form 
can  be  written  as 


f(s1,s2)  = 


a ~S1/91  ' s2/02 
-s—  e 


1 s 2 U 

c4-™) 


where  s^  > 0,  s2^  0,  and  Ig(z)  = £ 

k 0 


(z/2) 


is  the 


modified  Bessell  function  of  the  first  kind  and  order 
zero.  Here  a > 0,  d ^ 0,  and  a + d = 1. 

The  density  (4.10)  has  mean  vector 


Qj/a 


02/a 


C4.ll) 


covariance  matrix 


I - C°i j ) 


dy^z 


du1u2 


(4.12) 


and  correlation  p=d  with  0 p < 1.  The  marginal  dis- 
tributions of  Sj  and  S2  are  exponential  with  means  y^ 
and  y2  respectively.  A generalization  of  (4.10)  due  to 
Paulson  admits  of  correlation  values  -.25  <_  p < 1. 
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The  variates  (S ^ ,S2),  given  by  (4.10),  can  be 
readily  simulated.  Let  {X—  , j = l,2,...}  be  independent 
identically  distributed  exponential  with  mean  0^  for 
i=l,2  and  let  N denote  a geometric  random  variable  with 
density  function 


Pr  [N=i ] = p q,  i=l,2,. . 


where  p=d  and  q=a. 


Then 


csi,s2)  c xxj , x 2j) 


(4.13) 


has  the  bivariate  exponential  distribution  (4.10)  (see 
Downton  (1970)).  The  simulation  proceeds  by  simulating  N 
and  then  adding  that  number  of  i.i.d.  exponentials  accor- 
ding to  (4.13). 

It  would  be  easy  to  construct  situations  in  which 
the  service  times  are  negatively  correlated.  In  order 
that  we  may  readily  consider  this  case  we  require  a bi- 
variate exponential  distribution  which  is  conveniently 
simulated.  Such  a distribution,  which  includes  the  dis- 
tribution in  (4.10)  as  a special  case  is  due  to  Paulson 

® 1 

(1973).  The  actual  simulation  of  variates  S = Q from 

s2 

that  distribution  is  effected  through 


§*,  = + Vil2  + V1  V2  *3  + ; 


(4.14) 
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here  Xj  is  a 2-vector  of  independent  exponential  vari- 
ates with  mean  vector 


and  the  Vj  are  random  2x2 


0 o] 

1 o' 

0 oj  ’ 

L°  0. 

0, 


matrices  which  take  on  values  in  the  set  { 

[o  l]  ’ [2?]  > with  probabilities  a,b,c,d  respectively 

All  the  X j 1 s and  V^'s  are  mutually  independent.  Note 
that  eventually  the  product  HVj  will  result  in  a matrix 
of  zeros  and  so  with  probability  one  is  represented 
by  a finite  sum  (Kesten  (1973),  Kohberger  (1975)).  It 
turns  out  that  (4.13)  is  really  a special  case  of  (4.14). 

The  bivariate  random  variable  S in  (4.14)  has 
mean  vector 
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u = 
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V1 
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^2 
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(4.15) 


and  covariance  matrix 
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* 2 
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4.2.3  Simulation  Results  and  Interpretation 

Simulated  results  are  presented  in  this  section 
for  three  cases  of  interstage  storage  capacity:  (A) 

infinite  (q=°°)  , (B)  zero  (q=l)  , and  (C)  one  (q=2)  . For 
the  infinite  interstage  storage  case  results  will  be 
given  for  two  stages  in  series  for  various  values  of 
correlation  and  for  two  through  twenty-five  stages  in 
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series  for  correlation  equal  unity.  The  latter  depicts 
how  adding  stages  might  affect  system  performance  given 
correlation  p>0;  more  precisely,  it  provides  an  envelope 
within  which  system  performance  will  vary  since  for  a 
fixed  number  of  stages  and  utilization  correlation  unity 
provides  an  extremum  and  correlation  zero  provides  another. 
In  each  of  the  three  cases  we  allow  infinite  storage  be- 
fore stage  one. 

In  the  ordinary  case  in  which  the  correlation  be- 
tween paired  service  times  is  zero  a few  steady  state  re- 
sults are  available  for  comparison  purposes. 

We  have  taken  the  mean  arrival  rate  to  be  unity 
and  so  the  steady  state  utilization,  v,  at  stage  i is 
simply  the  mean  service  time  y^.  It  will  suffice  for  our 
purposes  to  take  y^  = y2  - y since  similar  steady  state 
behavior  will  obtain  for  y^  f y2 . Furthermore,  there  do 
not  seem  to  be  many  results  available  for  purposes  of  com- 
parison for  Cases  B and  C when  y^  j*  y2-  Tor  X = 1 , our 
system  performance  measure  of  mean  waiting  time  (queueing 
plus  service) , is  equivalent  to  the  expected  number  in 
the  system. 

Case  A.  k stage  queues  in  series,  infinite  interstage 
storage . 

(Graphs  are  labeled  k Q for  k-queues)  . 

Steady  state  results  for  k stages  in  series  with 
no  correlation  between  pairs  of  service  times  are  avail- 
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able  (Saaty ( 196 1) } and  we  have  that  the  expected  number 
of  customers  at  each  stage  is  v/(l-v)  and  kv/(l-v)  in 
the  system. 

The  following  diagram  provides  for  two  stages 
in  series  the  mean  waiting  time  at  the  second  stage, 
W<2).  for  v = 0.75  and  p = -0.25  , 0 , 0.50  , and  1.0.  In 
this  case  the  mean  waiting  time  at  the  first  stage  is 
independent  of  p since  no  blocking  occurs  and  hence  it 
suffices  to  examine  the  mean  waiting  time  at  the  second 
stage  to  determine  the  effects  of  correlated  service 
times.  In  some  of  the  simulation  results  to  follow  we 
replicate,  many  times,  runs  of  much  shorter  length;  here 
we  choose  to  illustrate  the  mean  waiting  time  as  a func- 
tion of  n with  one  very  long  run.  Long  runs,  such  as 
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this  one,  may  be  considered  as  being  composed  of  replica- 
tions of  smaller  runs  where  the  starting  condition  of  a 
new  replication  is  the  ending  value  of  the  previous  repli- 
cation (Conway(  1963)).  From  the  diagram  it  is  clear  that 
each  graph  tends  to  stabilize  for  increasing  n in  accor- 
dance with  the  law  of  large  numbers. 

The  numbers  adjacent  to  the  values  of  p in  the 
diagram  are  the  mean  waiting  times  for  the  second  stage 
and  mean  waiting  times  in  the  system  after  100,000  service 
completions.  We  point  out  that  for  p=0  the  mean  waiting 
times  at  stages  one  and  two  are  2.99  and  3.05  respectively 
and  these  are  in  close  agreement  with  the  expectation  of 
3.0  for  this  utilization.  For  p^0  we  see  here  a bonus 
attached  to  positive  correlation  in  service  times  since 
system  performance  improves  with  increasing  correlation. 

On  the  other  hand,  system  performance  deteriorates  with 
negative  correlation. 

Each  illustration  like  this  one  has  a starting 
condition  based  on  the  mean  waiting  time  from  a pilot  run 
and  then  we  omit  the  waiting  times  of  the  first  1,000  cus- 
tomers in  the  actual  computations  shown. 

Figure  3 gives  system  performance  for  different 
values  of  v and  for  p=0,l.  These  graphs  are  intended  to 
show  that  there  is  no  discemable  effect  due  to  correla- 
tion p>0  at  utilization  v = 0.6  but  as  v increases  from 
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0.6  to  0.9  a definite  trend  appears. 

Figure  4 shows  the  ratio  of  mean  time  in  the 
system  for  various  values  of  p to  the  mean  waiting  time 
in  the  system  at  p = 0 (Paulson  and  Beswick(  19 73)).  These 
kinds  of  graphs  are  based  on  an  average  of  100  replica- 
tions of  1,400  service  completions  (after  an  initializa- 
tion of  400  service  completions  were  discarded) . 

The  solid  lines  in  Figures  4,  7,  and  9 depict 
a smoothed  fit  to  the  actual  data.  Sampling  variation, 
of  course,  precludes  the  possibility  of  obtaining  such 
a smooth  fit  without  extremely  long  runs  or  extensive 
replications  but  each  curve  was  spot-checked  to  ascertain 
whether  or  not  the  fit  was  spurious.  In  no  case  was  any 
substantial  deviation  recorded. 

Now  we  show  how  these  effects  are  consistent 
with  Conolly's  (1968)  results  for  the  case  of  v=0.9  and 
correlation  of  p=1.0  between  the  service  times  in  the  two 
stages.  Conolly  showed  for  his  single  server  queueing 
system  where  the  ratio  of  service  time  to  the  interarrival 
time  was  constant  for  all  n,  that  for  a utilization  of  0.9 
(the  ratio)  the  mean  waiting  time  (queueing  plus  service) 
was  2.71.  For  service  time  independent  of  interarrival 
time  the  steady  state  expectation,  for  this  utilization, 
is  9.0.  The  interarrival  time  and  service  time  in  Conolly's 
system  are  perfectly  correlated  whereas  in  our  system  the 
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two  service  times  are  perfectly  correlated.  It  is  clear 
from  equation  (4.2)  that  the  correlation  between  the  in- 
terarrival time  at  the  second  stage  and  the  service  time 
there  is  less  than  one  and  so  the  improvement  in  system 
performance  for  our  system  should  be  less  than  Conolly's 
(an  elementary  derivation  shows  the  correlation  to  be  vp 
or  0.9  in  this  case).  We  see  from  Figure  3 that  the  mean 
waiting  time  at  the  second  stage  after  100,000  customers 
is  4.13  and  indeed  the  improvement  is  less. 

In  Rgure  5a  we  show  the  mean  waiting  time  as  a 
function  of  n for  five  stages  in  series  where  the  service 
times  are  equal  at  each  stage.  The  graphs  are  labeled 

W1-  1 corresponding  to  the  mean  waiting  time  at  stage  k, 

—(21 

k=l,2,...,5.  We  see  that  the  mean  waiting  time  W1-  1 , for 
the  second  stage,  is  consistent  with  the  results  in  Rgure 
3.  The  results  for  W^^,  W^^  , and  W^  suggest  that  fur- 
ther improvements  in  system  performance  occur  over  the  p=0 
case  but  the  effect  seems  to  approach  a limit.  The  number 

— f kl 

in  parenthesis  to  the  right  of  Wv  1 is  the  ratio 
k _n 

l Wl  J / (kv/ (1 -v) ) , 2 < k < 5.  Rgure  5b  shows  this  ratio 
i = 1 

for  two  through  twenty-five  stages  in  series  for  p=l.  These 
results  were  obtained  by  extending  recursive  formulae  (1) , 
(2),  and  (3).  In  this  extreme  case  of  correlation,  adding 
stages  has  an  effect  on  system  performance  which  depends 
markedly  on  the  utilization  rate;  e.g.,  for  v=0.7  system 
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performance  is  improved  through  the  first  four  stages 
and  then  is  reduced.  A utilization  of  0.9  gives  rise 
to  much  improved  system  performance  through  twenty* five 
stages . 

Cases  B and  C.  Two  stage  queues  in  series,  finite  (in- 
cluding zero)  interstage  storage. 

For  these  cases  the  utilization  is  effec- 
tively reduced  in  value  (Saaty (1961)).  The  maximum  effec- 
tive utilization  is  vmax  = (q+l)/(q+2)  where  the  queue  in 
stage  two  is  limited  to  a length  of  q-1  units.  We  consider 
the  cases  q=l  and  q=2. 

Figure  6 shows  the  mean  waiting  time  at 
the  first  stage  for  q=l  and  several  values  of  v.  For  this 
case  each  customer's  waiting  time  at  the  second  stage  is 
simply  the  service  time  there  so  we  are  concerned  only  with 
the  waiting  time  process  at  stage  one.  Figure  7 shows,  for 
stage  one,  the  ratio  of  mean  waiting  time  at  stage  one  with 
p^O  to  the  mean  waiting  time  at  stage  one  with  p=0. 

Steady  state  results  fcr  the  mean  number 
of  customers  in  the  system,  L,  for  p = 0,  q=l  and  with  util- 
ization v are  given  in  Morse  (1958);  we  have  that 

L = 4v(2-v2)/((2+v) (2  - 3v) ) . (4.17) 

Tor  v = 0.4,  0.5,  and  0.6  and  for  p = 0,  the 
observed  (expected)  values  of  L are  1.55  (1.53),  2.87  (2.80), 
and  7.80  (7.57)  respectively.  The  observed  values  are  from 
Figure  6. 
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For  p/0  , again  we  see  a dramatic  effect 
in  system  performance.  System  performance  deteriorates 
as  the  correlation  p increases  through  positive  values 
and  improves  as  p decreases  through  negative  values.  FFence 
when  there  is  no  storage  alloAved  before  stage  two  the  de- 
parture from  independence  results  in  significantly  differ- 
ent steady  state  behaviors,  especially  as  the  value  v 


is  approached. 


Finally,  we  consider  the  case  q=2.  Figure 


8 shows  the  mean  waiting  time  as  a function  of  n at  the 
first  stage  and  the  ending  value  for  the  mean  waiting 
time  in  the  second  stage.  The  mean  waiting  time  at  stage 
two  was  very  stable  for  all  values  of  n so  those  values 
will  not  be  illustrated.  The  effect  for  v=0.6  is  in  the 
same  direction  as  for  q=l  but  reverses  as  v increases  so 
that  for  values  close  to  v the  change  in  system  perfor- 
mance  is  consistent  with  the  q=°°  case;  that  is,  improve- 
ment for  p>0  and  deterioration  for  p<0.  Figure  9 shows 
the  effect  in  this  case  for  p/0. 

4.2.4  Spectral  Analysis  of  } and 

In  this  section  we  review  briefly  the  theory  of 

spectral  analysis,  show  the  sample  power  spectra  of  the 

time  series  {W^^}  and  } and  finally  apply  a non- 

n n 

parametric  test  to  the  ratio  of  certain  estimated  power 
spectra . 
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Several  authors  give  complete  accounts  of  the 
theory  and  application  of  spectral  analysis;  e.g.  see 
Anderson  (1971)  and  Jenkins  (1961).  Fishman  and  Kiviat's 
(1967)  paper  on  the  analysis  of  simulation  generated  time 
series  is  also  of  direct  interest.  Our  objective  here  is 
to  review  enough  of  the  fundamentals  of  spectral  analysis 
to  motivate  the  nonparametric  test  to  be  presented  later. 

In  general  we  take  (Zt,  teT)  to  be  a stochastic 
process  and  we  let  (zt,  teT)  denote  a sequence  of  obser- 
vations from  the  process;  the  sequence  is  referred  to  as 
a realization  of  the  process  or  simply  as  a time  series. 

We  take  the  index  set  T to  consist  of  discrete,  equispaced 
time  points.  From  the  time  series  we  seek  to  describe  the 
underlying  process.  Seldom  can  we  determine  the  form  of 
the  multivariate  distribution  which  generated  the  realiza- 
tion and  most  often  we  must  make  simplifying  assumptions 
even  to  describe  any  of  the  distribution's  moments. 

We  assume  that  the  process  is  in  a particular  state 
of  equilibrium  where  the  first  and  second  moments  are  inde- 
pendent of  time.  Therefore, 

E(ZJ  = y (4.18) 

l 

and 


c°vCZt,Zt.k)  - Yk 


(4.19) 
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Yg  = /g  f(cj)  dw.  (4.24) 

Thus  the  power  spectrum  can  be  considered  as  a decomposi- 
tion of  the  variance  at  different  frequencies. 

To  get  sample  results  which  are  statistically  con- 
sistent we  do  not  estimate  the  spectrum  at  a particular 
frequency  but  instead  estimate  the  average  power  about  the 
frequency  of  concern.  The  average  power  corresponds  to 
weighting  the  autocovariances  in  the  time  domain  and  we 
typically  estimate  f(w)  with  the  truncated  estimate 
~ i 5? 

f(Wj)  * -[Xg  Cg  ♦ 2 J ck  cos  2Tru>jk]  (4.25) 

where  uij  * j/(2m),  j-0, 1,2, . . . ,m  and  the  weights  Ak, 
k=0,l,2, . . . ,m,  form  a so-called  lag  window.  We  choose 
the  Blackman-Tukey  "hamming”  window, 

\k  ■ 0.54  ♦ 0.46  cos  ffk/m,  k*0 , 1 , 2 , . . . ,m.  (4.26) 

In  (4.25),  the  sample  autocovariances  cm+j>  cm+2»**,» 
are  omitted  since,  for  m sufficiently  large,  they  should 
contribute  little  information.  As  a result,  only  m auto- 
covariances need  be  calculated  and  savings  in  computation 
may  be  considerable.  Considerable  care  must  be  used  when 
selecting  m,  however,  because  too  large  a value  will  in- 
crease the  variance  of  the  estimates  and  too  small  a value 
will  not  give  enough  resolution. 
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Next  we  examine  several  sample  power  spectra  as- 
sociated with  the  simulated  waiting  times  for  the  two- 
server  infinite  interstage  storage  case.  We  take  the  sim- 
ulated values  (wCi))»  n*l , 2 , . . . ,N ; i*l,2,  to  be  time  series 
where,  as  before,  W^^  is  the  total  waiting  time,  queueing 
plus  service,  of  customer  n at  server  i.  Figure  10  shows 
a portion  of  the  sample  spectra  for  {W^1^}  and 
n«l,2, . . . ,2000,  and  for  correlation  values  of  p=0,  0.25, 
0.50,  and  1.0.  Utilization,  v,  is  0.90.  The  2000  sample 
values  were  chosen  from  the  end  of  a simulation  run  of 
length  30,000  to  ensure  that  any  possible  effects  of  start- 
up conditions  were  eliminated.  After  making  several  pilot 
runs,  m in  equation  (4.25)  was  set  equal  to  400.  R>r 
p*0.50  and  p=1.0  in  Figure  10  it  is  obvious  that  the  waiting 
times  at  the  second  server  give  rise  to  different  spectra 
than  the  waiting  times  at  the  first  server. 

Since  the  integral  of  the  power  spectrum  measures 
the  variance  of  the  process  and  the  area  under  the  sample 
spectrum  should  be  indicative  of  the  sample  variance,  we 
see  that  the  effect  of  positive  correlation  is  to  reduce 
the  variance  of  the  waiting  time  process.  Again  this  is 
consistent  with  Conolly's  results  (1968)  for  the  single 
server  system  in  which  a customer's  service  is  completely 
determined  by  the  length  of  the  interarrival  interval  sep- 
arating himself  and  his  predecessor.  R>r  a utilization 
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of  0.9,  Conolly's  system  reduces  the  steady  state  variance 
of  the  waiting  time  from  81,  for  the  classic  M^l/1  system, 
to  1.16;  the  sample  variance  associated  with  the  waiting 
times  at  server  2 in  figure  lOd  for  p*1.0  is  2.05  (the  sam- 
ple variance  associated  with  the  waiting  times  at  server 
1 is  60.1).  Recall  from  Section  4 that  the  condition  p=1.0 
for  the  correlation  between  a customer's  service  times  at 
the  two  servers  is  equivalent  to  a correlation  of  v,  or  0.9 
in  this  case,  between  his  interarrival  time  and  service 
time  at  the  second  server.  Therefore,  the  reduction  in 
variance  is  consistent  with  Conolly's  results  since  the 
corresponding  correlation  in  his  system  is  one. 

Next  we  develop  a nonparametric  test  for  the  hypo- 
thesis that  f ^ (co)  * f(2)  (w) , 0 a)  _<0.5,  where  f^(w) 
represents  the  power  spectrum  at  frequency  u>  associated 
with  the  time  series  } , n-l,2,...,N;  i»l,2.  The 

Blackman -Tukey  "hamming"  lag  window  in  (4.25)  gives  rise 
to  spectral  estimates  which  are  not  independent  and  so  we 
employ  the  notion  of  equivalent  independent  estimates 
(Jenkins (1961))  which  implies,  for  this  window,  that  esti- 
mates are  approximately  independent  if  they  are  about 
5/ (4m)  cycles  apart.  Since  the  estimates  in  (4.25)  are 
separated  by  a basic  frequency  of  l/(2m)  cycles,  this 
spacing  of  5/(4m)  cycles  amounts  to  taking,  as  independent, 
those  estimates  which  are  separated  by  an  interval  of  2.5 
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times  the  basic  frequency.  Since  the  waiting  times  are 
not  normally  distributed  and  the  assumption  of  normality 
is  implicit  in  the  development  of  equivalent  independent 
estimates  we  take  this  spacing  of  2.5  times  the  basic  fre- 
quency simply  to  be  a rough  guide.  Actually,  the  normality 
assumption  is  more  critical  for  making  distributional  as- 
sumptions about  the  spectral  estimates  than  for  the  usage 
here.  To  select  a practical  spacing  and  to  reduce  any 
possible  effects  of  the  normality  assumption,  we  take  es- 
timates at  the  frequencies  j/(2m),  j-1,4,7, . . . , to  be  ap- 
proximately independent  (the  spacing  here  is  3 times  the 
basic  frequency).  Therefore,  of  the  401  estimates  in  each 
spectrum  partially  illustrated  in  Figure  10,  we  take  134 
estimates  at  the  frequencies  j/800,  j-1,4,7, . . . ,399,  to  be 
approximately  independent. 

Now  for  each  approximately  independent  estimate  we 
can  regard  the  ratio  f ^ (u>)/f^  (w)  as  a Bernoulli  trial 
(greater  than  unity  ot  less  than  unity)  and  under  the  null 
hypothesis  of  homogeneity  of  the  two  spectra,  we  can  take 
as  a test  statistic  the  number  of  ratios  which  are  less 
than  unity.  Figure  11  shows  the  ratio  for  p*0  and  p-0.25. 
Of  the  134  approximately  independent  ratios  in  Figure  11a » 
64  are  less  than  unity  and  in  Figure  lib,  43  of  the  134 
ratios  are  less  than  unity.  Uider  the  null  hypothesis,  a 
ratio  greater  than  unity  is  as  equally  likely  as  a ratio 
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less  than  unity  and  taking  a normal  approximation  to  the 
implied  binomial  distribution,  we  have  a probability  of 
0.31  associated  with  observing  64  or  fewer  ratios  less  than 
unity  in  11a  and  a probability  of  .001  associated  with  ob- 
serving 43  or  fewer  ratios  less  than  unity  in  lib.  Although 
not  illustrated,  the  results  for  p=0.50  and  p=1.0  are  even 
more  conclusive:  for  p*0.50,  23  of  the  134  approximately 

independent  estimates  are  less  than  unity  and  for  p*1.0, 
none  of  the  ratios  are  less  than  unity.  Therefore,  for  the 
case  presented  here  we  have  good  statistical  evidence  that 
the  power  spectra  associated  with  the  waiting  times  at  each 
server  are  not  homogeneous  for  correlation  p>0.  We  expect 
similar  results  for  other  values  of  correlation,  utiliza- 
tion, and  interstage  storage  to  obtain. 

In  the  next  section  we  study  in  a similar  way  single 
server  queues  (not  in  tandem)  which  have  correlated  inter- 
arrival and  service  processes. 

4 . 3 Single  Server  Queues  with  Correlated  Interarrival  and 
Service  Processes 

In  this  section  we  are  concerned  with  a single  server, 
first-come,  first-served  queueing  system  where  we  depart  from 
the  usual  assumption  of  independence  by  taking  a customer's  in- 
terarrival interval  and  subsequent  service  times  to  be  correlated 
according  to  the  bivariate  exponential  distribution  described  in 
4.2.2.  Here  the  customer's  interarrival  interval  is  measured  be- 
tween his  arrival  time  and  that  of  his  predecessor.  It  is  assumed 
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that  customers  from  an  infinite  population  arrive  at  a single 
stage  according  to  a Poisson  process  with  rate  X,  which  we 
take,  without  loss  of  generality,  to  be  unity;  an  unlimited 
queue  is  allowed. 

As  in  the  previous  section,  the  system  performance  measure 
is  taken  to  be  the  mean  waiting  time  per  customer  and  in  this 
section  we  show  a formula  to  recursively  compute  the  waiting 
time  per  customer.  We  give  simulated  results  which  show  the 
effect  of  correlation  and  apply  the  nonparametric  test  of  the 
last  section  to  show  that  the  effect  is  indeed  statistically 
significant. 

Denote  by  Tn  the  time  between  arrival  epochs  of  customers 

cn_2  and  cn  to  the  queue  and  let  cn  experience  the  service  time 

Sn,  n-1,2,3,...  . The  sequence  of  interarrival  times  (Tn)  and 

service  times  (S „}  for  different  customers  are  both  assumed  to 
n 

be  independent;  for  customer  cn  we  assume  the  r.v.  (Tn»sn)  has 
the  bivariate  exponential  distribution  given  in  (4.10)  with 
p^-1  (to  be  associated  with  Tn)  and  U2”v  C^or  sn) » 0 < v < 1» 
so  that  the  steady  state  utilization  is 

R>r  customer  cn  we  define  wn  to  be  the  waiting  time,  ex- 
cluding service,  and  Wn  to  be  the  total  waiting  time,  n-1,2,...; 
the  following  diagram  illustrates  the  definitions,  it  is  ob- 


vious that  a recursive  formula  for  Vi 


is 


"n  * Vl  * Vl-  lf  Tn+i  « 


®n.l  » 


lf  Tn*l  i wn  . 
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Figure  13  shows  the  ratio  of  mean  time  in  the  system 
for  various  values  of  p to  expected  waiting  time  in  the  system 
at  p=0  and  we  see  that  the  effect  of  nonzero  correlation  is 
greatest  for  large  utilizations.  In  every  case  simulation  re- 
sults were  in  close  agreement  with  known  results  (for  p = 0 and 


Next  we  examine  sample  power  spectra  asoociated  with  the 
waiting  time  process  and  test  the  hypothesis  fg(w)  = fp(w), 

0 _5  uj  » where  fg(u>)  is  the  power  spectrum  associated  with 
the  waiting  time  process  at  p=0  and  likewise,  fp(u)  corresponds 
to  p^O. 

We  take  the  simulated  values  (Wn>,  n*l  ,2  , . . . ,2000  , to  be  a 
time  series  and  Figure  14  shows  a portion  of  the  sample  power 
spectra  for  p=0,  0.50,  and  1.0;  utilization  is  0.70.  The  spec- 
tra appear  different  and  we  suspect  that  the  variances  of  the 
waiting  time  processes  decrease  with  positive  correlation  due 
to  the  relation  of  the  illustrated  graphs.  In  fact,  the  corre- 
sponding simulated  waiting  time  series  have  variances  as  follows: 
p*0,  4.763;  p*0.50,  2.141;  and  p*1.0,  0.514.  (The  expected 

variance  for  p=0  is  5.444  C^orse  1958)  and  Conolly's  model  for 
p*l  gives  rise  to  a variance  of  0.613.)  Coupled  with  the  simu- 
lated result  that  p*-0.25  leads  to  a variance  of  10.451  we  see 
that  the  effect  of  positive  correlation  is  to  reduce  the  variance 
of  the  waiting  time  process  and  negative  correlation  causes  an 
increase . 
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As  a function  of  the  approximately  independent  estimates 

A A 

defined  in  4.2.4,  Figure  15  shows  the  ratios,  fg  (oo)  / fp  (go)  , for 
p-0.5  where  the  caret  signifies  that  we  are  using  sample  esti- 
mates. Recall  that  we  regard  each  ratio  as  a Bernoulli  trial 
(greater  than  unity  or  less  than  unity)  and  under  the  null  hy- 
pothesis of  homogeneity  of  the  two  spectra  we  take  as  a test 
statistic  the  number  of  ratios  which  are  less  than  unity.  In 
the  figure,  22  of  the  134  approximately  independent  ratios  are 
less  than  unity  which  is  very  strong  evidence  that  the  null  hy- 
pothesis is  false;  for  other  cases  of  p the  number  of  ratios 
less  than  unity  are:  p=-0.25,  87;  p=0.25,  43  and  p=1.0,  none. 

Additionally  we  applied  the  test  to  fQ  25 (w)/£q  50 (u)  and  got 
38  of  the  134  approximately  independent  ratios  less  than  unity. 
We  reject  the  implied  null  hypotheses  in  all  cases  and  conclude 
that  the  waiting  time  process  as  a function  of  p leads  to  dif- 
ferent power  spectra. 

R>r  interest's  sake  we  investigated  the  system  under  study 
for  a different  bivariate  exponential  distribution.  Primarily 
due  to  the  ease  with  which  the  variates  can  be  simulated,  we 
chose  the  bivariate  exponential  distribution  of  Marshall  and 
Olkin  (1967).  If  the  r.v.'s  U,  V and  W are  independent  exponen- 
tials with  parameters  Xj,  X2  and  X12,  respectively,  then  the 
bivariate  r.v.  (T,  S) , where  T=min(U,  W) , S-min  (V,  W)  , has  the 
indicated  distribution.  It  can  be  shown  that  T and  S have  means 
1/ (Xx  ♦ X^2)  and  1/(X2  + X^2)  , respectively,  and  the  correlation 
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between  the  two  is  + A-,  + A^)  • ^°r  t^ie  ProPer  choice 

of  parameters  then,  we  can  generate  a bivariate  r.v.  where 
E[T]=1,  E[S]  = v and  the  correlation  is  p. 

Figure  16  shows  comparative  results  for  the  Marshall  and 
Olkin  model  and  the  Wicksell-Kibble  special  case  in  (4.10). 
Although  not  illustrated,  we  applied  the  homogeneity  test  to 
the  two  sample  spectra  shown  with  the  result  that  65  of  134 
approximately  independent  ratios  (fw_K(aO _0(<jo)  ) were  less 
than  unity.  The  probability  is  approximately  0.37  of  observing 
65  or  fewer  ratios  less  than  unity  if  the  hypothesis  fw  k((d)  = 
is  true  so  we  fail  to  reject  it.  We  expect  similar  re- 
sults for  other  values  of  correlation  and  utilization  to  obtain. 

4 . 4 Summary 

We  have  extended  the  work  of  Paulson  and  Beswick  (1973) 
for  the  effect  of  dependent  exponential  service  times  on  the 
system  performance  of  tandem  queueing  systems.  We  assumed  that 
each  queue  has  a single  server  and  the  service  discipline  is 
first-come  first-served. 

For  the  two  stage  queueing  system  we  derived  recursive 
formulae  for  the  total  waiting  time  per  customer  at  each  queue 
for  the  cases  of  zero,  one,  and  infinite  interstage  storage. 
Simulation  results  for  the  mean  waiting  time,  under  the  assump- 
tion of  correlated  service  times,  showed  that  the  system's  be- 
havior is  quite  sensitive  to  departures  from  the  traditional 


assumption  of  mutually  independent  service  times,  especially 
at  higher  utilization  rates.  For  the  case  of  infinite  inter- 
stage storage,  mean  waiting  time  is  reduced  by  positive  corre- 
lation and  increased  by  negative  correlation.  This  change  is 
reversed,  however,  for  zero  interstage  storage  and  depends  on 
the  value  of  the  utilization  rate  for  the  case  where  interstage 
storage  equals  unity.  By  using  spectral  analysis  and  a non- 
parametric  test  applied  to  the  sample  power  spectra  associated 
with  the  simulated  waiting  times  we  showed  that  the  effect  is 
statistically  significant;  in  addition  we  showed  that  the  vari- 
ance of  the  waiting  time  process  is  reduced  for  positive  corre- 
lation. 

We  showed  in  a precise  way  how  the  two  stage  queueing  sys- 
tem with  dependent  service  times  and  infinite  interstage  storage 
is  related  to  a single  server  system  with  interdependent  arrival 
and  service  processes;  an  interpretation  by  Conolly  (1968)  for  a 
special  type  of  this  latter  interdependence  is  shown  to  be  use- 
ful in  suggesting  why  the  mean  and  variance  of  the  waiting  time 
process  are  affected  by  correlated  service  times. 

Par  correlation  equal  unity  and  infinite  interstage  storage, 
results  were  shown  for  two  through  twenty- five  stages  in  series; 
these  results  provide  an  envelope  within  which  system  performance 
will  vary  since  for  a fixed  number  of  stages  and  utilization 
correlation  unity  provides  one  extremum  and  correlation  zero 
provides  another. 
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Additionally,  we  studied  single  server,  single  stage 
queueing  systems  wherein  a customer's  interarrival  interval 
and  subsequent  service  times  are  governed  by  a bivariate  ex- 
ponential distribution.  For  the  case  of  unlimited  storage 
capacity  we  showed  that  positive  correlation  leads  to  reduced 
mean  waiting  times  and  negative  correlation  increases  the  mean 
waiting  times,  both  more  so  at  higher  utilizations.  The  results 
were  shown  to  be  statistically  significant.  We  also  showed  that 
the  variance  of  the  waiting  time  process  is  reduced  for  positive 
correlation  and  increased  for  negative  correlation.  Briefly  we 
investigated  the  effect  of  using  the  Marshall  and  Olkin  (1967) 
bivariate  exponential  distribution;  results  were  similar. 


I 
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PART  5 


DISCUSSION  AND  CONCLUSIONS 

Here  we  have  investigated  various  bivariate  non-normal 
distributions  and  showed  several  areas  to  which  they  could  be 
applied.  Of  principal  concern  were  bivariate  negative  binomial 
and  gamma  distributions.  In  Part  2 we  discussed  numerous  theor- 
etical aspects  and  fitted  the  distributions  to  data.  Part  3 il- 
lustrated some  bivariate  approaches  for  analyzing  aircraft  op- 
erations and  maintenance  data  and  lastly,  Part  4 showed  how 
bivariate  exponential  distributions  could  be  applied  to  queueing 
systems  with  certain  kinds  of  correlation.  Summary  results  were 
given  in  each  part;  next,  we  highlight  these  results  and  suggest 
areas  for  future  research. 

One  bnb  distribution  was  obtained  by  convolving  the  Paulson- 
Uppuluri  (1972)  bivariate  geometric  distribution  which  is  de- 
fined by  a certain  characteristic-functional  equation.  The  dis- 
tribution has  six  parameters  and  admits  of  positive  or  negative 
correlation  and  linear  or  nonlinear  regression  functions.  Shown 
were  the  moments  to  order  two  and  for  special  cases,  the  regres- 
sion function,  a recursive  formula  for  the  cell  probabilities, 
a method  of  moments  parameter  estimation  technique,  the  likeli- 
hood equations,  the  differential-difference  equations  and  for 
maximum  likelihood  estimation,  a necessary  relationship  for  the 
parameters.  Certain  analogous  properties  were  shown  for  a dual 
bivariate  gamma  distribution.  For  both  of  these  distributions 
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areas  for  future  research  include  determining  an  infinitely 
divisible  representation  for  the  full  model  (all  six  para- 
meters) and  developing  parameter  estimation  techniques.  Both 
areas  represent  sizeable  tasks. 

We  discussed  another  bnb  distribution  which  resulted  by 
reducing  a certain  trivariate  negative  binomial  distribution 
with  independent  marginals.  Moments  to  order  two  and  the  re- 
gression function  were  provided.  This  distribution  has  mar- 
ginals whose  characteristic  functions  are  of  the  form 
. -v . 

[l+Sfl-e11)]  , whereas  the  aforementioned  distribution's 

marginal  characteristic  functions  are  like  [1  + 0^  (l-e1*)  ]"v,  i*l  ,2 . 
Although  not  investigated  here  we  suspect  that  bnb  distribu- 
tions which  allow  0 and  v to  vary  would  have  practical  value. 

That  these  bnb  distributions  should  be  useful  was  illus- 
trated by  analyzing  sample  data  sets,  some  with  negative  corre- 
lation and  nonlinear  regression. 

Another  siajor  effort  centered  on  forming  bivariate  r.v.  's 
related  to  particular  aircraft  operations  and  maintenance  prob- 
lems. For  a random  sample  we  showed  that  the  negative  binomial 
distribution  could  be  used  to  adequately  describe  demands  for 
aircraft  spare  parts  for  single  time  periods  (univariate)  and 
adjacent  time  periods  (bivariate).  An  application  for  fly-away 
kits  was  discussed.  Additionally,  we  investigated  aircraft 
abort  data  for  single  six  month  periods  and  adjacent  six  month 
periods  and  showed  how  the  univariate  and  bivariate  negative 


99 


binomial  distributions  fitted  the  data.  The  P values  assoc- 
2 

iated  with  the  x goodness-of- fit  test  ranged  between  0.01 
and  0.85  with  the  average  being  about  0.36.  With  these  abort 
data  we  illustrated  a possible  way  to  examine  the  effect  of  a 
noted  event  on  an  item's  performance.  Here  we  were  interested 
in  the  effect  of  overhaul  on  an  aircraft's  performance  when 
measured  by  aborts  for  six  month  periods.  Basically  the  method 
involved  comparing  two  bivariate  distributions,  one  defined  for 
r.v. 's  on  either  side  of  the  event  (overhaul)  and  the  other  de- 
fined for  a similar  r.v.  not  separated  by  the  event.  We  used 
the  regression  functions  to  compare  the  sample  distributions. 

For  our  data,  aborts  increased  after  overhaul  but  since  our 
analysis  was  limited  in  certain  ways  we  were  unable  to  conclude 
that  the  rise  was  due  to  overhaul.  We  intend  to  investigate 
this  application  more  in  the  future. 

The  last  part  showed  the  effect  of  certain  correlated  r.v.' 
on  the  system  performance  of  tandem  and  single  stage  queueing 
systems.  A bivariate  exponential  distribution  was  used.  In 
both  cases  we  assumed  that  arrivals  were  according  to  a Poisson 
process,  the  service  discipline  was  first-come,  first-served 
and  a single  server  was  available. 

For  the  two  stage  tandem  queueing  system  we  showed,  via 
simulation,  that  the  mean  waiting  time  is  quite  sensitive  to 
departures  from  the  traditional  assumption  of  mutually  indepen- 
dent service  times,  especially  at  higher  utilizations.  For 
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the  case  of  infinite  interstage  storage,  mean  waiting  time 
is  reduced  by  positive  correlation  and  increased  by  negative 
correlation.  This  change  is  reversed,  however,  for  zero  inter- 
stage storage  and  depends  on  the  value  of  the  utilization  rate 
for  the  case  where  interstage  storage  equals  unity.  A result 
by  Conolly  (1968)  was  shown  to  be  useful  in  explaining  the 
effect  for  the  infinite  interstage  storage  case. 

For  single  stage  queueing  systems  where  a customer's 
interarrival  interval  and  subsequent  service  times  are  corre- 
lated we  showed  that  positive  correlation  reduces  mean  waiting 
time  and  negative  correlation  increases  mean  waiting  time. 

The  storage  area  was  assumed  to  be  infinite. 

By  using  spectral  analysis  and  a nonparametric  test  ap- 
plied to  the  sample  power  spectra  associated  with  certain  sim- 
ulated waiting  times  we  showed  the  effect,  in  both  cases,  to 
be  statistically  significant. 


* 


PART  6 


REFERENCES 


Abramowitz,  M.  and  Stegun,  I. A.,  eds.  (1964).  Handbook  of 
Mathematical  Functions.  Superintendent  of  documents , 

U.S.  Government  Printing  Office,  Washington,  D.C. 

Anderson,  T.  W.  (1971).  The  Statistical  Analysis  of  Time 
Series , John  Wiley  6 Sons , Inc.  , New  York. 

Arbous , A.  G.  and  Kerrich,  J.  E.  (1951).  "Accident  Statistics 
and  the  Concept  of  Accident-Proneness, " Biometrics,  Vol.  7, 
340-432. 

Arbous,  A.  G.  and  Sichel,  H.  S.  (1954).  "New  Techniques  for 
the  Analysis  of  Absenteeism  Data,"  Biometrics,  Vol.  41, 
77-90. 


Arnold,  B.  C.  (1967).  "A  Note  on  Multivariate  Distributions 
with  Specified  Marginals,"  J.  Amer.  Statist.  Ass.,  Vol. 

62,  1460-1. 

Bates,  G.  E.  and  Nevman,  J.  (1952).  "Contributions  to  the 

Theory  of  Accident  Proneness,  I,"  University  of  California 
Publications  in  Statistics,  215-54. 

Bhat,  U.  N.  (1969).  "Queueing  Systems  with  First-Order  De- 
pendence," Opsearch,  Vol.  6,  No.  1,  1-24. 

Block,  H.  W.  (1975).  "Physical  Models  Leading  to  Multivariate 
Exponential  and  Negative  Binomial  Distributions*"  Rensselaer 
Polytechnic  Institute,  0.  R.  and  S.  Research  Paper  No. 
37-75-P3. 

Boswell,  M.  T.  and  Patil,  G.  P.  (1970).  "Chance  Mechanisms 
Generating  the  Negative  Binomial  Distributions,"  published 
in  Random  Counts  in  Models  and  Structures , Vol.  1,  edited 
by  G.  P.  Patil,  Penn.  State  Univ.  Press,  3-22. 

Burke,  P.  J.  (1956).  "The  Output  of  a Queueing  System,"  Opera- 
tions Research,  Vol.  4,  699-704. 

Clark,  J.  R.  (1972).  Properties  of  a New  Multivariate  Geometric 
and  Negative  Binomial  Distribution  with  Possible  Applications, 
M.3.  tnesis,  Univ.  of  Tenn. 

Conolly,  B.  W.  (1968).  "The  Waiting  Time  Process  for  a Certain 
Correlated  Queue,"  Operations  Research,  Vol.  16,  1006-1015. 


101 


102 


Conolly,  B.  W.  and  Hadidi,  N.  (1969).  MA  Correlated  Queue," 

J.  Appl . Prob.,  Vol . 6,  No.  1,  122-136. 

Conolly,  B.  W.  and  Hadidi,  N.  (1974).  "A  Comparison  of  the 
Operational  Features  of  Conventional  Queues  with  a Self- 
Regulating  System,"  Appl,  Statist. , Vol.  18,  41-53. 

Conway,  R.  W.  (1963).  "Some  Tactical  Problems  in  Digital 
Simulation,"  Management  Science,  Vol.  16,  No.  1,  47-61. 

Cross,  K.  (1970).  A Gradient  Projection  Method  for  Constrained 
Optimization,  Report  No.  tf-1746.  Oak  Ridge  National  Labora- 
tory, Oak  Ridge,  Tennessee. 

Dade,  M.  (1973).  Examples  of  Aircraft  Scheduled-Maintenance 
Analysis  Problems,  The  Rand  Corp. . R-l299-pR. 

Downton,  F.  (1970).  "Bivariate  Exponential  Distributions  in 
Reliability  Theory."  J.  Royal  Statist.  Society,  Series  B, 

Vol.  32,  408-417. 

Draper,  N.  R.  and  Smith,  H.  (1966).  Applied  Regression  Analysis, 
John  Wiley  6 Sons,  Inc.,  New  York. 

Edwards,  C.  B.  and  Gurland,  J.  (1961).  "A  Class  of  Distribu- 
tions Applicable  to  Accidents,"  JASA,  Vol.  56,  503-517. 

Faucett,  W.  M.  and  Gilbert,  R.  D.  (1966).  Characteristics 
of  Demand  Distributions  for  Aircraft  Spare  Parts,  Research 
and  Engineering  Departments  (ERR-FW-512) , General  Dynamics, 
Fort  Worth  Division. 

Fishman,  G.  S.  and  Kiviat,  P.  J.  (1967).  "The  Analysis  of 
Simulation-Generated  Time  Series,"  Management  Science, 

Vol.  13,  525-557.  

Gibbons,  J.  D.  (1971).  Nonparametric  Statistical  Inference, 
McGraw-Hill  Book  Company,  New  York. 

Guldberg,  A.  (1934).  "On  Discontinuous  Frequency  Functions 
of  Two  Variables,"  Skand.  Aktuar. , Vol.  17,  89-117. 

Harris,  R.  (1968).  "Reliability  Applications  of  a Bivariate 
Exponential  Distribution."  Operations  Research,  Vol.  16, 

No.  1,  18-27. 

Hawkes,  A.  G.  (1972).  "A  Bivariate  Exponential  with  Applica- 
tions to  Reliability,"  J. Royal  Statist.  Society,  Series  B, 
Vol.  34,  129-131. 


Holgate,  P.  (1964).  "Estimation  for  the  Bivariate  Poisson 
Distribution,"  Biometrika,  Vol.  51,  241-S. 

Jackson,  R.  R.  P.  (1954).  "Queueing  Systems  with  Phase  Type 
Service,"  Operations  Research,  Vol.  5,  109-120. 

Jenkins,  G.  M.  (1961).  "General  Considerations  in  the  Analysis 
cf  Spectra,"  Technometrics , Vol.  3,  No.  2,  133-166. 

John,  F.I.  (1963).  "Single  Server  Queues  with  Dependent  Ser- 
vice and  Interarrival  Intervals ,"  J.  Soc.  Indust.  Appl. 
Math.  , Vol.  11,  526-534.  — 

Johnson,  N.  L.  and  Kotz,  S.  (1969).  Discrete  Distributions, 
Houghton  Mifflin  Co.,  Boston. 

Jury,  E.  I.  (1964).  Theory  and  Application  of  z-Transform 
Method,  John  Wiley  § Sons,  Inc.,  New  York. 

Kemp,  C.  D.  (1970).  "Accident  Proneness  and  Discrete  Distri- 
bution Theory,"  published  in  Random  Counts  in  Scientific 
Work , Vol.  2,  edited  by  G.  P.  Patil,  Penn.  State  Univ. 

Press , 41-66. 

Kendall,  M.  G.  and  Stuart,  A.  (1973).  The  Advanced  Theory  of 
Statistics , Vol.  2,  Hafner  Publishing  Co.,  New  York. 

Kesten,  H.  (1973).  "Random  Difference  Equations  and  Renewal 
Theory  for  Products  of  Random  Matrices  ,"  Acta  Mathematica, 
Vol.  131,  207-248. 

Kibble,  W.  F.  (1941).  "A  Two-Variate  Gamma  Type  Distribution," 
Sankhya,  Vol.  5,  137-150. 

Kohberger,  R.  C.  (1975).  "On  Certain  Multivariate  Exponential 
Distributions ,"  Ph.D.  Thesis,  Rensselaer  Polytechnic  Insti- 
tute, Troy,  New  York. 

Krishnamoorthy , A.  S.  and  Parthasarty,  M.  (1951).  "A  Multi- 
variate Gamma  Type  Distribution  ,"  Ann.  Math.  Statist., 

Vol.  22,  549-557.  

Lindley,  D.  V.  (1952).  "The  Theory  of  Queues  with  a Single 
Server  ."  Proc.  Camb.  Phil.  Soc.,  Vol.  48,  277-289. 

Lundberg,  0.  (1940).  On  Random  Processes  and  their  Application 
to  Sickness  and  Accident  Statistics,  Uppsala:  Almquist  and 

Wiksell. 


J 

I 


I 


104 


Mann,  N.  R. , Schafer,  R.E.,  and  Singpurwalla , N.D.  (1974). 

Methods  for  Statistical  Analysis  of  Reliability  and  Failure 
Data,  John  Wiley  5 Sons,  Inc. .New  York 

Mardia,  K.  V.  (1970).  Families  of  Bivariate  Distributions, 
Hafner  Publishing  Co. , New  York. 

Marshall,  A.  W.  and  Olkin,  I.  (1967).  "A  Multivariate  Expon- 
ential Distribution."  J.Amer.  Statist.  Assoc.,  Vol.  62, 

30-44. 

Morse,  P.  M.  (1958).  Queues,  Inventories  and  Maintenance, 

John  Wiley  § Sons , Inc. , New  York. 

Paulson,  A.S.  (1973).  "A  Characterization  of  the  Exponential 
Distribution  and  a Bivariate  Exponential  Distribution," 
Sankhya,  Series  A,  Vol.  35,  69-78. 

Paulson,  A.  S.  and  Beswick,  C.  A.  (1973).  "The  Effect  of  De- 
pendent Exponential  Service  Times  on  Queues  in  Series," 
Department  of  Operations  Research  and  Statistics,  Rensselaer 
Polytechnic  Institute  Research  Report  No.  37-73-P2. 

Paulson,  A.  S.  and  Uppuluri,  V.  R.  (1972).  "A  Characterization 
of  the  Geometric  Distribution  and  a Bivariate  Geometric 
Distribution,"  Sankhya,  Series  A,  Vol.  34,  297-300. 

Saaty,  T.  L.  (1961).  Elements  of  Queueing  Theory  with  Appli- 
cations , McGraw-Hill,  New  York. 

Subrahmaniam,  K.  (1966).  "A  Test  for  "Intrinsic  Correlation" 
in  the  Theory  of  Accident  Proneness,"  J.  Royal  Statist. 
Society,  Series  B,  Vol.  28,  180-189. 

Subrahmaniam,  ffocherlakata  and  Subrahmaniam,  Kathleen  (1973). 

"On  the  Estimation  of  the  Parameters  in  the  Bivariate 
Negative  Binomial  Distribution."  J.  Royal  Statist.  Society, 
Series  B,  Vol.  35,  131-146. 

Titchmarsh,  E.  C.  (1964).  The  Theory  of  Functions,  London, 
Oxford  University  Press. 

Whittaker,  E.  T.  and  Watson,  G.  N.  (1965).  A Course  of  Modern 
Analysis , Cambridge  University  Press. 

Wicksell,  S.  D.  (1933).  "On  Correlation  Functions  of  Type  III," 
Biometrika,  Vol.  25,  121-133. 

Youngs,  J.W.T.,  Geisler,  M.A.  and  Brown,  B.B.  (1955).  The 
Prediction  of  Demand  for  Aircraft  Spare  Parts  Using  the 
Method  of  Conditional  Probabilities,  The  Rand  Corp.,  RM- 1413. 


'■  •.  • - - - - —L  . 

y 


( 


APPENDIX  A 

A NEW  MULTIVARIATE  NEGATIVE  BINOMIAL  DISTRIBUTION 

In  this  appendix  we  show  Clark's  (1972)  derivation  of  a 
new  multivariate  negative  binomial  distribution.  The  distri- 
bution results  by  convolving  a multivariate  geometric  distri- 
bution. 

The  multivariate  analogue  of  (2.13)  is 

4>(T)  - <KT)E[<KTV)  ] (A.  1) 

where 

T = (tj  » *2  * * * * »*n^  * 

n Pi  Ui  -l 
*(T)  * n [l+TlI_(l-e  J)]  1 

j«l  Fj 

and  V is  the  set  of  all  n dimensional  diagonal  matrices  of 
zeroes  and  ones  (2n  matrices  in  all).  Here  veV  assumes  a par- 
ticular value  with  probability  a , £ a ■ 1 and  a + aT  < 1, 

v i 0,  I where  £ is  the  zero  matrix  and  I is  the  identity  matrix. 
Equation  (A.l),  which  can  be  written  as 

♦ CT)  - *(T)  l *v4>(Tv),  (A. 2) 

veV  - ~ 

defines  a multivariate  geometric  distribution. 

It  can  be  shown  that 

it.  , 

♦(0,...,0,tm,0 0)  - [l+0m(l-e  m)]"1  (A. 3) 

where 
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and  v*eV  such  that  mth  diagonal  term  is  one  (2n  * matrices 
in  all).  Therefore  the  marginals  are  geometric  and  it  fol 
lows  that 


- V 

Var  (X  1 = 6 (1  + 0 ) , m=l,2,...,n. 
v nr  m v nr  ’ ’ ’ ’ 

In  (A.  2)  if  ay  = 0 for  v t (),  I we  have 
♦ (T)  * MT)(a+d<J>(T)) 


(A. 4) 


(A. 5) 


where  a is  associated  with  ag  and  d with  a^.  Convoluting  as  in 
(2.25)  yields 

♦V(T)  ' Ia*(T)]v[l-d*(T)]‘v 

- [a^(T)]V[l+(^)d<i;(T)  + (V21)(dlKT))2+...]  (A. 6) 

and  inverting  gives  the  multivariate  negative  binomial  dis- 
tribution 


where 


gv(xi,x2,...,xn)  - C[nFn.1(v+x1,v+x2, 

n 

v+x  v ;d  n (1  -p  - ) ) J , 

n i*=l  1 


n 

n 


V+X.-l  X. 

( 3 )(l-pi)v(pi)  3 

xj  3 3 


f 
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. x 


and  _F 


nFn-l(al an ’ bl bn-l;z) 


, . ? (ai>r"tan)i  *: 

1 ik  U>i>  -Oj  T 


(A.  7) 


Clark  concluded  by  showing  that 


E[Xm>  - a0» 


and 


Var  [Xm]  = v9mCl+0m).  m=l,2 n; 


he  also  showed  that  the  distribution  is  infinitely  divisible 
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APPENDIX  B 


RECURSIVE  FORMULA  FOR  WAITING  TIMES 

In  this  appendix  we  show  how  a set  of  recursive  formulae 
for  the  waiting  times  can  be  constructed  for  any  number  of 
queues  in  series  where  interstage  storage  is  unlimited. 

Referring  to  the  following  diagram  for  cn  and  cn+1's 
queueing  and  service  times  at  stages  two  and  three  (a  continua- 
tion of  the  diagram  preceding  equation  (4.1)  in  the  text),  we 
see  that  cn+^’s  interarrival  time  at  stage  three  is 
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Similar  to  equation  (4.3)  in  the  text,  cn+j's  waiting 
time  at  the  third  stage  is 


W 


n+1 


(3)_ 


= < 


n+1 , 3* 


if  Tn.l,3iW„C3) 


W t3)-T  _+S 

n n+1, 3 n+1, 3’ 


if  T ,,  ,<W 
n+1 , 3 n 


(3) 


Comparing  Tn  + 1 2 anc*  Tn  + 1 3 we  ^ave,  in  general,  for 


cn+1's  interarrival  time  at  stage  i,  i=2,3,..., 
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Similarly,  comparing  wn  + 1^1'*»  wn+l^^  and  Wn+1^^  gives 
a general  recursive  formula  for  cn+1's  waiting  time  at  stage 
i,  i-1,2, . . . , 
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Thus,  we  can  obtain  the  recursive  formulae  for  any  number 
of  queues  in  series. 
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TABLE  1.  BIVARIATE  DATA  SETS 
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Probability  of  exceeding  computed  x value  with  indicated  degrees  of  freedom 
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Guldberg-Bates -Neyman  bivariate  negative  binomial  distribution  of  (2.6) 
estimates  are:  6=2.8956,  0*1.5854.  There  results  y2=17.0  and  P=0.20  f 


TABLE.  3.  BWB(aJo,o)pJ^y)  MODEL. 
EXPECTED  CELL  FREQUENCIES  FOR 
ARBOUS-S'.CREL  DATA  (2h8  WORKERS) 
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Clark's  bivariate  negative  binomial  distribution  of  (2.26).  ML  estimates  are: 
a=0 . 0370 , p=0.1010,  q=0.0970,  0=1.5480. 
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TABIC  4.  G-B-N(«t,e,w  ) MODEL  . 
observed  and  expected  cell  FREQUENCIES 
FOR  BATES-NEVMAN  DATA  (1286  WORKERS). 
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OBSERVEO  AND  EXPECTED  CELL  FREQUENCIES 
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Paulson-Uppuluri  bivariate  geometric  distribution 

of  (2.36).  ML  estimates  are:  a«0.1338,  b=0.0191, 

c»0.0330,  p*0.1860,  q-0.4516.  There  results 
202.4  and  Pi0  for  df»87. 
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TAftL£  C>  8M8  la^b.O.p.q,  v ] MODEL. 

OBSERVED  ANO  EXPECTED  CELL  FREQUENCIES 
FOR  BATES-NETMAM  DATA  (123  6 WORKERS) . 
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TABLE.  T BNBCa.o.o,  p,q,  y)  MOOEL1  . 
OBSERVED  ANO  EXPECTED  CELL  FREQUENCIES 
FOR  BATES-NETMAM  DATA  (1286  WORKERS). 
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^Clark’s  bivariate  negative  binomial  distribution  of 

(2.26).  ML  estimates  are:  a-0.3725,  p-0.3086, 

0.6299,  v*1.1654.  There  results  x^*195.4  and  P«0 
for  df-89. 
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TABLtl  0.  BNB-TR(e,w'.iw’i,L'3)  MODEL. 
OBSERVED  AND  EXPECTED  CEIL  FREQUENCIES 
FOR  BATES-NEWN  DATA  (1286  WORKERS). 
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sults x^"152.6  and  P-0  for  df-82. 
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TABLE  v observed  and  expecteo  ceu.  FREQUENCIES 
OF  DEMAND  FOR  72  AIRCRAFT  PARTS 
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TABLE  U.  MODEL*. 

OBSERVED  AND  EXPECTED  CELL  FREQUENCIES 
OF  FLIGHT  ABORTS  FOR  lOT  AIRCRAFT. 
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^aulson-Uppuluri  bivariate  geometric 
distribution  of  (2.36).  ML  estimates 

A ^ A A 


are:  a=0, 

q-0 . 3299 . 
P*0 . 12  for 


b=0 . 6820 , c-0. 3179,  p=0.1655 
2 

There  results  x *10.2  and 
df*6. 
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TOTAL  ABORTS  TOR  PERIOD  1,A, 


TABLE.  13.  BNBCa.c^o.p^.y)  MODEL  . EXPECTED  CELL  FREQUENCIES 
OFTOTAL  ABORTS  FOR  20  3 AIRCRAFT. 
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TABLE  1^+  8NB  (a ,o,o,  p,  w)  MODEL1. 
OBSERVEO  AND  EXPECTED  CELL  FREQUENCIES 
OF  TOTAL  ABORTS  FOR  203  AIRCRAFT. 

NO  \ntervening  overhaul. 
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TOTAL  ABORTS  FOR  PERIOD  I^A, 


LClark's  bivariate  negative  binomial  distribution 
of  (2.26).  ML  estimates  are:  a=0.6795,  p=0.6827, 

q»0.6971,  C“1 • 7701.  There  results  x^=46.1  and  P= 
0.65  for  df=50. 


TABLE  15  OBSERVED  CELL  FREQUENCIES  OF  TOTAL  ABORTS 
FOR  38 T AVRCRAFT.  INTERVENING  OVERHAUL. 
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TABLE  17  BNB (a,o,o,p.q.w>)  MODEL1. 
OBSERVED  AMD  EXPECTED  CELL  FREQUENCIES 
OF  TOTAL  ABORTS  FOR  J07  AIRCRAFT. 
INTERVENING  overhaul. 


^lark ' s,  bivariate  negative  binomial  distribution 
of  (2.26).  ML  estimates  are:  a=0.7470,  p=0.7526, 

q*0.7939,  v*1.7359.  There  results  x^*10?.0  and 
P-0.12  for  df»91. 
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CORRELATION 

FIGURE  4.  RATIO  OF  MEAN  WAITING  TIMES,  2Q.q 
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FIGURE  10.  A PORTION  OF  WAITING  TIME  SPECTRA  ,2Q,i/=0  9, 


FIGURE  12.  MEAN  WAITING  TIMES,  v=0.7. 
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